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operations._2) 
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strength.  Though  much  still  remains  to  be  learned,  it  is  now  possible  to  make  qualitative 
engineering  predictions  on  bearing  strength  as  a result  of  the  comprehensive  laboratory  and 
field  research  program  directed  at  defining  the  tensile  strength,  compressive  strength,  flexural 
strength,  elastic  modulus,  creep  behavior,  and  their  relationship  to  temperature  and  salinity. 
Bearing  strength  analysis,  using  the  CF.t.-developed  finite-element  computer  codes,  can  now 
be  based  on  the  material  properties  corresponding  to  the  actual  temperature  and  salinitv 
gradient  existing  across  the  vertical  profile  of  the  ice  plate. 
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METRIC  CONVERSION  CHART 


1 i 
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°F 

Celsius  (C) 

(°F-32)  X 0.556 

INTRODUCTION 

!n  recent  years  operational  use  ol 
ice-covered  waters  in  polar  regions  has 
greatly  increased  to  coincide  with  the 
advent  of  a world-wide  energy  crisis.  No 
longer  can  remote  polar  locations  be  con- 
sidered too  inhospitable  or  uneconomical 
to  search  for  or  harvest  energy  re- 
sources. In  particular,  Naval  Antarctic 
operations  continue  to  depend  on  sea-ice 
airfields  and  roads  to  provide  heavy 
cargo  and  logistics  support  for  United 
States  activities.  Therefore,  a better 
understanding  of  the  behavior  of  sea-ice 
sheets  under  aircraft  and  vehicular 
loading  is  required  to  ensure  continued 
safe  operations.  To  this  end  the  Civil 
Engineering  Laboratory  (CEL)  was  initi- 
ally engaged  in  performing  experiments 
and  developing  analyses  in  which  sea-ice 
sheets  were  assumed  to  behave  elastic- 
ally  (1,  2,  3).  The  culmination  of  this 
elastic-properties  research  on  sea  ice 
was  the  development  of  an  elastic  finite- 
element  computer  code  that  has  the 
capability  of  modeling  temperature- 
dependent  properties  as  a function  of  ice 
sheet  thickness  (*♦). 

However,  it  is  known  that  sea  ice 
behaves  elastically  only  for  small  to 
moderate  load  magnitudes  applied  over 
short  time  durations.  Consequently, 
parked  aircraft,  cargo  storage,  over-ice 
platforms  or  stations,  and  even  snow 
overburden  cannot  be  correctly  analyzed 
by  elastic-behavior  solutions  to  deter- 
mine their  actual  effect  on  ice  sheets. 
To  define  the  range  of  greater  load 
magnitudes  and  longer  time  durations, 
CEL  conducted  another  series  of 
laboratory  and  field  experiments  on  the 


creep  behavior  of  sea  ice.  The  theory  of 
linear  viscoelasticity  was  coupled  with 
the  finite-element  technique  to  provide 
more  realistic  characterization  of  sea- 
ice  sheet  behavior  (5). 

The  objectives  of  this  report  are:  (1) 
to  summarize  the  CEL  elastic  and 
viscoelastic  experimental  programs;  (2) 
to  review  both  finite-element  techniques 
developed  by  CEL  for  ice-sheet  analysis; 
and  (3)  to  present  bearing  capacity  limits 
for  both  short-  and  long-term  ice-sheet 
operations. 


SEA-ICE  PHYSICS 

Before  proceeding  with  mechanical 
properties  experimental  programs  on  sea 
ice,  it  is  necessary  to  review  physical 
material  properties  of  sea  ice  to  under- 
stand its  structure  and  p.  )vide  insight 
into  its  behavior  under  load.  Sea  ice  is 
sensitive  to  several  influential  param- 
eters, e.g.,  salinity,  temperature,  crystal 
orientation,  density,  and  seasonal 
changes.  A typical  phase  diagram  is 
presented  in  Figure  1 that  shows  the 
relationship  between  ice,  brine,  and  solid 
salts  as  a function  of  temperature  (6). 
For  any  specified  temperature  this 
diagram  gives  the  relative  volume  of 
brine,  ice,  and  solid  salt  in  the  system,  as 
well  as  depicting  the  temperatures  of 
crystallization  of  the  major  solid  salts 
found  in  sea  ice.  This  figure  assumes  that 
ion  distribution  is  the  same  as  in 
standard  seawater;  thus,  any  selective 
brine  drainage  would  cause  changes  in 
this  standard  sea-ice  phase  diagram. 
Once  these  relative  volumes  of  major 
sea-ice  components  are  established,  it  is 
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important  to  understand  how  sea  ice 
forms,  both  its  macrostructure  and 
microstructure. 

Macrostnicture  of  Sea  ke 

A sea-ice  sheet  begins  to  grow  as 
minute  spheres  that  quickly  change  into 
thin,  wafer-like  discs,  and  then  into  star- 
like  crystals  which  float  on  the  surface. 
These  stars  grow  rapidly,  overlap,  and 
start  to  freeze  together  into  a thin  ice 
skim.  If  the  water  should  remain  calm, 
this  ice  skim  would  have  a completely 
vertical  c-axis  orientation.  However, 
since  the  water  surface  can  never  be 
perfectly  calm,  a number  of  star  crystals 


get  tipped  and  caught  with  an  inclined  c- 
axis. 

At  this  stage,  the  growth  rate  is 
determined  by  the  thermal  gradient  in 
the  young  ice  sheet.  Those  ice  crystals 
with  a vertical  c-axis  grow  slower  than 
the  surrounding  horizontal  c-axis  crys- 
tals. Thus,  any  tendency  toward  aniso- 
tropic growth  produces  geometric 
selection,  finally  eliminating  all  but 
those  crystals  having  a horizontal  c-axis. 
This  transition  takes  place  in  a very  thin 
surface  layer,  commonly  only  I to  2 in. 
thick.  Because  the  transition  layer  is 
thin,  it  often  sublimes  away  or  recrystal- 
lizes into  the  snow  cover.  Therefore, 
over  a single  ice-growth  season,  a sea- 


Figure  2.  Horizontal  thin  section  of  natural  seawater  ice  enlarged  to  show 
grain  and  subgrain  structure.  The  grid  is  1 cm  on  a side. 


ice  sheet  can  be  considered  to  be 
columnar  throughout  its  thickness. 
Columnar  sea  ice  means  that  it  has 
strong  crystal  elongation  parallel  to  the 
direction  of  heat  flow,  a definite  crystal 
orientation,  and  an  increase  in  crystal 
size  at  the  growth  front  (bottom)  of  the 
'ce  sheet.  Average  crystal  diameter  is 
aboui  O.ff  in.  at  a level  12  in.  below  the 
ice-sheet  surface. 

Microstnxrture  of  Sea  ke 

When  seawater  freezes,  the  dis- 
solved salts  are  rejected  by  the  ice, 
causing  a solute  concentration  in 
"pockets"  of  liquid  brine.  Each  sea-ice 
crystal  forms  as  a "packet"  of  plates 
separated  by  an  array  of  brine  pockets. 
As  shown  in  Figure  2,  a single  crystal  can 
be  identified  easily  under  crossed 


polaroids,  since  it  forms  a common 
extinction  unit.  There  is  a random 
orientation  of  crystal  c-axes  in  a hori- 
zontal plane,  but  within  each  crystal  the 
c-axis  has  only  one  horizontal  direction. 

Brine  pockets  arc  generally  tubular 
cells  standing  vertically  in  the  icc  sheet. 
The  distance  between  adjacent  rows  of 
b."ino  pockets,  parallel  to  the  c-axis,  is 
called  the  plate  spacing.  The  spacing 
does  vary  with  position  in  the  ice  sheet 
due  to  changes  in  growth  velocity,  but 
the  distance  is  on  the  order  of  0.02  in. 
On  the  other  hand,  the  distance  between 
individual  brine  pockets  at  a single  plate 
interface  varies  much  more  than  the 
plate  spacing,  primarily  due  to  tempera- 
ture. As  ice  temperatures  become 
colder,  sea  ice  will  have  smaller  and 
more  widely  spaced  brine  pockets  due  to 
freezing  of  the  water  and  precipitation 
of  solid  salts  from  the  brine. 


When  analyzed  structurally,  a 
potential  intragranular  failure  between 
cold  ice  plates  would  occur  on  a plane 
having  a larger  percentage  of  ice.  Thus, 
sea  ice  should  have  a much  higher 
strength  at  colder  temperatures.  A 
structural  theory  (7)  has  been  developed 
relating  sea-ice  strength  to  a material 
property  called  brine  volume,  u,  a func- 
tion of  ice  temperature  and  salinity. 
Consequently,  the  following  section  des- 
cribes the  different  ranges  of  test 
parameters  used  in  the  sea-ice  mechan- 
ical properties  experimental  program  and 
the  rationale  behind  their  selection. 

Laboratory  Sea-Ice 
Characteristics 

As  a naturally  occurring  material, 
sea  ice  :s  senailive  to  several 
environmental  parameters  (e.g.,  salinity, 
temperature,  seasonal  changes,  crystal 
orientation).  To  isolate  and  control 
these  parameters,  CEL  has  conducted 
comprehensive  mechanical  property 
studies  in  a cold-chambcr  environment, 
supplemented  by  appropriate  field 
experiments  to  provide  data  checkpoints. 
An  insulated  freezing  tank,  64  in.  by  101 
in.,  provided  unidirectional  ice  growth 
from  the  surface  to  duplicate  the  natural 
environment  for  sea-ice  growth  (8).  Since 
the  freezing  tank  was  not  equipped  for 
circulating  the  water  to  maintain  con- 
stant salinity,  the  ice  was  grown  from 
water  that  increased  in  salinity  as  brine 
was  rejected  during  ice  growth.  Such  a 
freezing  system  has  previously  been 
defined  as  a closed  system  (9). 

During  the  freezing  cycle  for  all  ice 
production  the  temperature  was  main- 
tained at  -29C.  Prior  to  freezing,  the 
water  was  prechilled  to  a temperature  of 
0.5C  to  1.5C.  During  the  prechill,  the 
water  was  occasionally  mixed  to  reduce 
chill  time  and  temperature  stratifica- 
tion. The  freezing  cycle  was  maintained 


until  there  was  an  ice  growth  of  22  to  24 
in.;  from  this  point,  the  temperature  of 
the  chamber  was  adjusted  to  the 
scheduled  test  temperatures.  Ice  that 
formed  after  the  beginning  of  the  test 
phase  was  excluded  from  the  sampling. 

The  salinity  of  the  water  under  the 
ice  sheet  at  the  end  of  the  freezing  cycle 
was  in  the  range  of  60  ppt  for  the  natural 
seawater  ice  and  9.5  ppt  for  the  brackish 
water  ice.  Freezing  and  testing  of  the 
ice  were  performed  in  the  same  cold 
chamber,  with  the  temperature 
fluctuating  +2.20  from  a given  test  or 
freezing  temperature.  The  temperature 
cycle  was  sufficiently  short  to  produce 
no  observable  effec*  on  the  steady-state 
temperature  of  the  test  specimens. 

Test  Temperature.  The  ice  test  tem- 
peratures, -4,  -10,  -20,  and  -27C,  were 
initially  chosen  to  bracket  the  precipita- 
tion points  of  the  two  major  salts  in 
normal  seawater:  sodium  sulfate  (Na2 
SO  4 • IOH2O)  with  a precipitation  point 
of  -8.2c,  and  sodium  chloride 
(NaCl  • 2H2O)  with  a precipitation 
point  of  -23. 3C  (see  Figure  1).  By 
performing  strength  property  tests  at 
temperatures  on  both  sides  of  the  major 
salt  precipitation,  the  effects,  if  any,  of 
solid  salt  reinforcement  could  be 
observed. 

Salinity.  Natural  seawater  was  used 
to  provide  a laboratory  sea  ice  of  7 to  9 
ppt  salinity  for  tensile,  flexural,  com- 
pression, and  shear  strength  experiments. 
Also,  tensile,  flexural,  and  shear  strength 
experiments  were  run  on  l-to-2-ppt 
brackish  ice,  a mixture  of  distilled  water 
and  natural  seawater.  Representative 
salinity  profiles  for  tanks  containing  both 
types  of  ice  are  shown  in  Figures  3 and  4. 
The  most  distinguishing  feature 
contrasting  the  vertical  profiles  of  the 
two  ice  types,  aside  from  difference  in 
actual  salinity  value,  is  a more  uniform 
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salinity  dsitribution  for  the  brackish  ice 
and  its  very  low  surface  salinity.  The  8- 
to-9-ppt  salinity  of  the  natural  seawater 
ice,  which  began  at  about  the  ^i-in. 
depth,  IS  fairly  typical  of  salinities 
rep.  rted  for  young  natural  sea  ice  (10). 
The  representative  horizontal  salinity 
profiles  in  Figure  ^ demonstrate  a fairly 
uniform  distribution,  and  they  are  in 
good  general  agreement  with  the  vertical 
profile  data  at  the  same  depth. 

Brine  Migration.  To  investigate 
effects  of  season  changes  and  detect 
brine  migration,  salinity  profiles  associ- 
ated with  thermal  activity  in  several 
laboratory-grown  sea-ice  blocks  were 
plotted  (Figure  5).  The  predominant  shift 
in  salinity  distribution  was  in  the  top  5 
in.  of  the  ice  block.  In  blocks  B and  C, 
the  decrease  was  8 ppt  and  5 ppt, 
respectively,  in  the  surface  stratum 
between  the  first  and  second  samplings. 
No  shift  was  observed  between  the 
second  and  third  samplings.  Prior  to  the 
second  sampling,  both  ice  blocks  had 
been  subjected  to  a condition  of  thermal 
stability;  that  is,  the  temperature  envi- 
ronment had  remained  constant. 
Between  the  second  and  third  samplings, 
bo. h ice  blocks  had  been  subjected  to  a 
condition  of  thermal  activity  by  the 
temperature  change.  The  brine  loss  in 
the  upper  stratum  is  weakly  indicated  by 
the  peaking  trends  picked  up  in  the  lower 
stratum  for  the  :>ccond  and  third  sam- 
plings; however,  no  appreciable  degree  of 
brine  migration  developed,  even  under 
thermal  activity,  during  the  800-to- 
.. 000-hour  observation  period.  The  early 
loss  of  salinity  in  the  first  few  inches  of 
the  surface  stratum  occurred  during  a 
condition  of  thermal  stability. 


0.920  to  0.925  gm/cm3,  while  that  for 
brackish  ice  was  found  to  be  slightly  less, 
A values  ranging  from  about 

0.900  to  0.915  grn/cm  A Representative 
density  profiles  are  presented  in  Figure 
6.  Both  liquid  displacement  and  dimen- 
sion measurement  techniques  were  used 
to  calculate  volumes. 

Grain  Size.  The  laboratory-grown 
sea  ice  was  very  similar  in  petrographic 
structure  to  natural  sea  ice,  having 
almost  identical  plate  spacing  (between 
U.3  and  0.5  mm).  However,  this  distance 
increased  to  almost  1.0  mm  for  brackish 
ice.  A pronounced  increase  in  grain  size 
with  depth  was  noted  for  natural 
seawater  ice,  while  grain  size  variation 
was  less  evident  for  brackish  ice  (see 
Figure  7). 

MECHANICAL  PROPERTIES  OF 
SEA  ICE 


Density.  The  jce  density  for  ruitural 
seawater  was  found  to  center  around 


In  formulating  a comprehensive 
experimental  program  for  measuring 
mechanical  properties  of  sea  ice,  every 
attempt  was  made  to  isolate  variable 
parameters  and  to  carefully  record 
testing  conditions.  In  addition  to  the  sea- 
ice  characteristics  described  previously 
It  is  known  that  sea  ice  under  load 
behaves  differe.  ly,  depending  on 
w ether  the  stress  is  applied  pcrpendicu- 
ar  or  parallel  to  crystal  growth. 
Laboratory  tests  were  conducted  to 
determine  the  strength  of  ice  specimens 
extracted  both  horizontally  and 
vertically  from  the  ice  block.  All  of  the 
strength  values  obtained  for  each  loading 
condition  are  related  to  the  square  root 
of  brine  volume,  which  is  expressed  as  a 
function  of  salinity  and  ice  temperature. 
The  following  sections  discuss  several 
strength  property  test  series  performed 
on  sea  ice. 
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Tensile  Strength 

Tensile  strength  is  best  determined 
by  direct  uniaxial  tension  testing.  Any 
substitute  for  the  direct  test  method, 
such  as  ring  tensile,  brazil,  or  flexural 
tests,  induce  complicated  stress  states 
within  the  sample  that  must  be  taken 
into  account  when  calculating  the  test 
results.  However,  uniaxial  tensile  testing 
for  brittle  materials  is  extremely  sensi- 
tive to  any  slight  stress  concentrations 
that  can  cause  premature  failure  and 
give  erroneous  results.  Therefore,  special 
care  was-  taken  in  designing  the  testing 
method  and  preparing  the  ice  specimens. 

Preparation  of  Dumbbell  Ice  Speci- 
mens. Tensile  tests  were  performed  on 
dumbbell-shaped  ice  specimens  to  reduce 
end-effect  errors.  First,  rough  cylindri- 
cal cores  were  extracted  from  sawed 
blocks  cut  out  of  laboratory-grown  sea 
ice  and  cut  into  right  circular  samples  by 
using  a cut-off  saw.  A special  alignment 
tray  to  hold  the  samples  during  cutting 
was  necessary  to  ensure  plane  and 
parallel  ends.  These  cylindrical  speci- 
mens were  cut  approximately  three  core 
diameters  long  to  provide  a test  section 
length  of  approximately  one  core  dia- 
meter. The  test  section  diameter  was 
one-half  of  the  core  diameter.The  transi- 
tion region  between  the  test  section  and 
enlarged  ends  must  avoid  stress  concen- 
trations due  to  short-radius  fillets.  For 
this  reason  a transition  radius  was  chosen 
almost  equal  to  the  core  diameter  - large 
enough  to  minimize  possible  stress 
concentrations.  The  tensile  ice 
specimens  were  machined  on  a regular 
metal-turning  lathe,  using  a single  com- 
pound cutting  tool  to  form  the  radius  and 
contact  seat  for  mating  with  the  test 
head.  The  entire  machining  setup  is 
shown  in  Figure  8.  After  machining,  the 
dumbbell  specimens  were  cold-soaked  at 
the  test  temperature  for  3 to  *»  hours 
prior  to  testing. 


Tensile  Pulling  System.  Most  past 
tensile  testing  procedures  have  favored 
freezing  the  dumbbell  ice  specimen  into 
metal  end  caps  (11).  As  an  alternative, 
CEL  developed  a tensile  pulling  system 
that  consisted  of  split  aluminum  end  cups 
with  a load-applying  pin  insert  and  a 
keeper  ring;  the  end  cups  were  machined 
to  receive  the  finished  ice  specimen. 
Both  the  gripping  heads  and  a dumbbell 
ice  specimen  are  shown  in  Figure  9.  The 
truncated  conical  section  of  the  ice 
specimen  bears  on  its  mating  cup  holder 
surface  and  transmits  the  uniaxial  load 
without  freezing  the  specimen  in.  Align- 
ment of  the  testing  machine  was  period- 
ically checked,  and  a universal  joint  was 
inserted  to  ensure  concentric  load 
application.  This  system  has  certain 
flexibility  to  readjust  the  specimen  after 
preloading  to  check  opposing  strain 
output.  In  most  cases,  however,  the 
specimen  slid  into  proper  position  guided 
by  the  beveled  bearing  surface.  This 
gripping  technique  made  testing  easier  to 
perform,  thus  allowing  a greater  number 
of  specimens  to  be  tested. 

To  check  the  alignment  of  the  tool, 
a model  ice  specimen  was  machined  out 
of  acrylic  plastic,  a photoelastic 
material.  The  fringe  pattern  symmetry 
shown  in  Figure  10  indicates  uniform 
tensile  stress  in  the  test  section. 
Photoelastic  stress  distributions  were 
checked  at  four  different  angles  around 
the  specimen.  Since  birefringence  is  a 
function  of  thickness,  a plot  of  the 
number  of  fringes  over  the  diameter 
should  describe  a semicircle,  another 
clicck  on  stress  uniformity.  In  addition 
to  this  photoelastic  test,  the  same 
acrylic  specimen  was  strain-gaged  at 
three  points,  120  deg  apart,  around  the 
circumference  of  the  test  section. 
Surface  strains  were  checked  under 
different  tensile  loads  for  uniformity. 
The  strain  measurements  from  the  three 
gages  were  found  to  be  within  3 to  5%  of 
each  other  for  three  separate  loading  and 
unloading  cycles. 
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Figure  8.  MKhining  tensile  specimen  with  compound  cutting  tool. 


Figured  Tensile  specimen  end  gripping  heeds. 


Figure  10.  Symmetric  fringe  pattern  in 
model  of  dumbbell  specimen  under 
tension. 


of  sea  ice  to  be  virtually  independent  of 
the  applied  stress  rate  within  the  above 
ranges  of  stress  loading  rate  (2).  Every 
sea-ice  specimen  loaded  in  tension  failed 
by  cleavage,  producing,  in  general,  a 
smooth  failure  plane  perpendicular  to  the 
axis  of  the  applied  load. 


Testing  Procedure.  All  tensile  ice 
specimens  has  a nominal  i.6-in.  diameter 
to  ensure  failure  without  end  effects; 
strength  calculations  were  based  on 
actual  diameter.  The  test  frame  (a 
10,000-lb  universal  testing  machine)  was 
located  in  the  cold  chamber  to  maintain 
a constant  specimen  temperature  during 
the  test.  A constant  crosshead  velocity 
of  0.5  in./min  was  maintained  throughout 
the  tests  to  provide  an  elastic  ice  re- 
sponse. This  produced  a stress  loading 
rate  of  15  to  35  psi/sec  for  horizontal 
specimens  and  25  to  <i5  psi/sec  for 
vertical  specimens.  A special  test  series 
indicated  the  horizontal  tensile  strength 


Laboratory  Results.  The  laboratory 
series  of  uniaxial  tensile  tests  on  sea  ice 
included  four  test  temperatures  (-4,  -10, 
-20,  and  -27C),  two  sample  orientations 
(horizontal  and  vertical),  and  two  salinity 
ranges  (7  to  9 ppt  and  1 to  2 ppt).  The 
sample  population  ranged  from  38  to  138 
for  each  of  the  16  possible  combinations 
of  the  above  parameters.  A summary  of 
the  tensile  strength  results  is  shown  in 
Figure  1 1 as  a function  of  the  sqaure 
roof  of  the  brine  volume,  v . The  average 
of  each  sample  group  is  plotted  with  a 
95%  confidence  level.  These  statistical 
limits  indicate  that  if  another  sample 
group  were  to  be  tested  with  the  same 
variables  held  constant,  there  would  be  a 
95%  chance  the  mean  of  that  sampling 
would  lie  within  the  confidence  interval. 
A curve-fitting  technique,  employing 
linear  regression  by  the  method  of  least 
squares,  is  used  to  determine  straight- 
line  equations  for  both  horizontal  and 
vertical  tensile  strengths,  o th  and  o.v, 
respectively,  where  o ^ is  found  in  Ib/in.  ^ 
and  V in  parts  per  thousand. 


0,,^  =118-  10.0 


(I) 


tv 


221*  - 12.7 


(2) 


By  the  least-squares  method,  different 
sample  groups  from  both  salinity  levels 
are  averaged  and  weighted  by  the 
amount  they  differ  from  the  mean.  Both 
Equations  1 and  2 are  plotted  on  Figure 
11.  The  y-intercept  of  these  equations 
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Figure  11.  Tensile  strength  versus  brine  volume. 


indicates  the  tensile  strength  of  ice 
having  zero  brine  content,  but  the  same 
structure  as  sea  ice.  The  negative  slope 
indicates  the  strength  of  ice  deteriorates 
as  V increases;  i.e.,  as  the  salinity 
increases  and/or  ice  temperature 
becomes  warmer.  Cc.mpicte  test  .esuits 
on  CEL's  tensile  strength  of  sea  ice  are 
presented  in  References  I and  2. 


salinity  level  (7  to  9 ppt).  The  sample 
population  ranged  from  46  to  81  for  each 
of  these  eight  combinations  of  sensitive 
parameters.  Loadhead  speed  was  kept 
constant  at  0.5  in. /min,  producing  a 
strain  rate  on  the  order  of  1 0*3  in./m. 
sec-‘.  Both  the  load  and  longitudinal 
deflection  were  monitored  during  each 
test. 


Compressive  Strength 

Laboratory  and  Field  Test  Results. 

Both  extensive  laboratory  and  Antarctic 
field  studies  were  performed  on  right- 
circular  cylinders  of  sea  ice  loaded 
elastically  in  compression  to  failure  (12). 
Each  specimen  had  a length-to-diameter 
ratio  of  2:1.  The  laboratory  test  series 
included  four  test  temperatures  (-4,  -10, 
-20,  and  -27C),  two  specimen  orienta- 
tions (horizontal  and  vertical),  and  one 


I he  field  test  series  included  just 
two  ice  temperatures  (-4  and  -IOC)  and 
only  the  vertical  crystal  orientation.  The 
salinity  range  (4  to  6 ppt)  was  consider- 
ably lower  than  normal  young  sea  ice, 
since  the  test  specimens  were  extracted 
from  year-old  sea  ice. 

While  every  attempt  was  made  to 
eliminate  salinity  as  a variable,  the 
average  salinity  varied  I to  2 ppt  within 
the  7-to-9-ppt  range.  In  Figure  12  both 
horizontal  and  vertical  compressive 
strength  are  plotted  as  a function  of  the 
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Figure  12.  Compressive  strength  versus  brine  volume. 


square  root  of  brine  volume,  v , again  having  no  biine  content  but  the  same  as 
using  the  linear  regression  technique  sea  ice. 

called  the  least-squares  method.  By  this  It  is  interesting  to  compare  the 
method  different  data  groups,  shown  in  least-squares  analyses  of  the  uniaxial 
Figure  12  as  individual  points  with  95%  tensile  strength  tests  with  those  of  the 
confidence  limits,  are  avcidgcu,  uniaxiai,  unconlined  compression 

weighted  by  the  amount  they  differ  from  experiments.  Two  significant  findings 
the  mean,  and  approximated  by  straight  can  be  determined  from  the  ratios  that 
lines  governed  by  the  following  equa-  show  the  relationships  between  (1) 

tions;  tensile  and  compressive  strengths  for 

each  specimen  orientation  and  (2) 
°ch  ■ 825  - 60.iyvr  (3)  horizontal  and  vertical  specimen  orienta- 

tion for  each  uniaxial  loading  mode. 

0 - 1 630  - 79  2 AT  iU)  *^*'^^*»  horizontal  and  vertical 

cv  ’ ' * cornpressive-to-tensile  strength  ratios 

are  approximately  7:1,  using  the  zero 
where  o ch^'’^  ^cv  horizontal  and  brine  values  for  comparison.  Second, 

vertical  compressive  strength  (lb/in.2)  of  both  the  compressive  and  tensile 

sea  ice,  respectively.  Again  the  y-inter-  vertical-to-horizontal  strength  ratios  are 
cept  of  the  above  equations  simply  approximately  2:1,  again  comparing  zero 
indicates  the  compressive  strength  of  ice  brine  volume  values. 
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Figure  13.  Horizontal,  compression  specimen  tested  at 
■20C.  showing  typical  the<ir  failure. 


Compressive  Failure  Modes.  The 
modes  of  failure  of  compressive  sea-ice 
specimens  vary,  depending  on  ice  temp- 
erature, specimen  crystal  orientation, 
end  re»itdirus,  and  strain  rate  during 
testing.  The  strain  rate  was  held  con- 
stant at  10 in./in./sec’*  for  these 
particular  tests.  Of  course,  compression 
itself  does  not  tend  to  form  or  propagate 
cracks  and,  thus,  cause  failure.  But  the 
same  does  not  hold  true  for  the  secon- 
dary tensile  and  shear  stresses  generated 
during  uniaxial  compression  (13).  In 
brittle  materials,  such  as  very  cold  sea 
ice,  tested  at  relatively  high  strain  rates, 
shear  fractures  inclined  to  the  compres- 
sion axis  are  quite  common  (Figure  13). 

Another  typical  failure  mode  of 
cylindrical  specimens  is  fracturing  along 
a plane  parallel  to  the  compression  axis. 
This  is  termed  axial  cleavage  and  is 


caused  by  secondary  tensile  stresses  in 
the  end  thirds  of  the  specimen.  Theore- 
tically, these  tensile  stresses  can  attain 
values  up  to  one-half  of  the  axial 
compressive  stress.  However,  they 
generally  are  suppressed  by  friction 
between  the  heads  of  the  testing 
machine  and  the  specimen.  This  end 
restraint  causes  an  abnormal  increase  in 
the  required  compressive  load  in  order  to 
build  up  the  tensile  stresses  sufficiently 
to  cause  failure.  End  confinement  can  be 
reduced  by  placing  low  friction  material, 
such  as  Teflon,  between  the  loadheads 
and  specimens.  Alternatively,  platen- 
specimen  matching  can  provide  the 
condition  of  equal  radial  strain  by 
forcing  the  quotient  of  Young's  modulus 
divided  by  Poisson's  ratio  to  have  the 
sane  value  for  both  specimen  and  platen 
materials  (11).  To  study  the  effects  of 


friction-reducing  materials  on  compres- 
sive strength,  CEL  undertook  a small- 
scale  laboratory  effort  discussed  in  a 
later  section. 

Almost  all  of  the  horizontal  ice 
specimens  tested  in  compression  failed  in 
shear,  with  higher  temperature  {-k  and 
-IOC)  specimens  failing  less  demonstrably 
than  those  at  lower  temperatures;  that 
is,  they  developed  a series  of  small, 
criss-crossing  shear  cracks,  but  remained 
intact.  The  greater  density  of  brine 
inclusions  at  these  higher  ice  tempera- 
tures may  create  more  points  of  stress 
concentration  which  ultimately  reduce 
the  load-carrying  capacity  of  the  speci- 
men without  its  destruction.  At  lower 
temperatures  a more  definitive  shear 
plane  develops,  because  the  sea  ice  has 
become  a more  homogeneous  crystalline 
material  after  increasing  salt  precipita- 
tion. 

Occasionally,  vertical  ice  specimens 
demonstrated  a cleavage  failure  due  to 
secondary  tensile  stresses  developing 
perpendicular  to  the  crystal  growth. 
These  tensile  stresses  tend  to  separate 
crystals  at  their  boundaries  or  within 
crystals  at  platelet  boundaries.  This 
failure  mode  may  govern,  particularly 
under  higher  load  rates  ( ~ 300  to  400 
psi/sec),  since  tensile  stresses  develop 
spontaneously  as  axial  compression  is 
rapidly  applied. 

Strain  Rate  Effects.  It  is  known  that 
different  loading  rates  alfcrt  compres- 
sion testing  significantly;  therefore,  a 
series  of  strain-rate  compressive 
strength  experiments  was  performed  on 
laboratory-grown  sea  ice.  The  tests  were 
conducted  at  two  different  ice  tempera- 
tures (-10  and  -20C)  on  right-cylindrical 
samples  cored  both  parallel  to  (vertical) 
and  perpendicular  to  (horizontal)  the 
direction  of  crystal  growth.  Average 
values  for  the  square  root  of  the  brine 
volume,  V , were  6.10  ppt  for  -IOC  and 


5.20  ppt  for  -20C.  The  cylindrical  ice 
cores  had  either  a 2-  or  3-in.  diameter 
and  were  cut  having  a length-to- 
diameter  ratio  of  2:1.  No  special 
friction-reducing  material  was  employed 
between  the  loading  platens  and  ice 
specimens.  Crosshead  travel  speeds 
ranged  from  a high  of  20  in. /min  to  a low 
of  0.05  in. /min.  In  Figures  14  and  15  both 
horizontal  and  vertical  compressive 
strengths  are  plotted,  respectively,  for 
two  ice  temperatures  and  one  complete 
range  of  loading  strain  rates.  All  of  the 
curves  show  a maximum  compressive 
strength  occurring  at  or  nc  .r  8.5  10"^ 

in./in.  sec  ■*,  which  corresponds  to  0.5 
in./min  crosshead  speed  cf  loading 
machine.  Shear  failure  planes  were  the 
predominant  feature  in  most  of  the 
horizontal  ice  specimens  (Figure  16); 
however,  at  the  faster  rates,  several 
specimens  failed  abruptly  by  cracking 
vertically  due  to  secondary  tensile 
stresses  (Figure  17).  This  latter  behavior 
pattern  was  ut:iversal  in  the  vertical  ice 
cores,  but  instead  of  a single  separation, 
they  fractured  at  grain  boundaries  to 
form  many  splintered  bundles  of 
columnar  ice  crystals. 

Failure  planes  develop  more  slowly 
corresponding  to  slower  loading  rates. 
Fvpn  though  the  ice  specimen  remains 
intact  at  slower  rates,  obviously  greater 
longitudinal  deformation  and  more 
microcracking  have  occurred  due  to 
plastic  flow.  Plastic  flow  is  the  sole 
cause  of  failure  below  3.4  x 10"^  in. /in. 
sec"*,  while  above  8.5  x 10"^  in. /in. 
’lec"*  sea-ice  specimens  behave  in  a 
brittle  manner.  As  the  load  is  rapidly 
applied,  a bursting  process  follows  the 
occurrence  of  brittle  cleavage  cracks 
and  the  linkage  of  grain  boundary  micro- 
cracks. There  appears  to  be  a transition 
phase  between  brittle  and  ductile 
response  at  loading  rates  ranging  from 
3.4  X 10^*^  to  8.5  X 10*^  in./in.  sec"*, 
where  both  slippage  and  brittle  cleavage 
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Figure  15.  Vertical  corrrpressive  strength  versus  strain  rate. 
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Figure  16.  Horizontal  compressive  ice  specimen 
tested  at  •20C  artd  loaded  at  a rate  of 
1 in./min  (1,7  x 10'^  in./in,  sec'M. 

cracks  may  occur  together.  Generally, 
microcracks  are  initiated  at  grain 
boundaries  and  are  followed  closciy  by 
abrupt  cleavage  cracking.  But  even  after 
cleavage  cracks  foim,  the  specimen  is 
able  to  sustain  additional  loading  and 
usually  attains  its  maximum  load  without 
explosive  rupture. 

Friction-Reduction  Tests.  Secondary 
tensile  stresses,  which  normally  would 


Figure  17.  Horizontal  comprenivt  ice  specimen 
tested  at  •20C  and  loaded  at  a rate  of 
5 in./min  (8.5  x 10'^  in. /in.  sec’^. 

cause  failure  at  lower  compressive 
strengths,  are  suppressed  due  to  friction 
building  up  between  the  ice  specimen  and 
the  loading  heads.  To  reduce  this 
friction,  some  type  of  material  is 
generally  placed  between  the  specimen 
and  platens  to  permit  lateral  movement. 
Four  different  materials  - cardboard 
(0.023-in.  thick),  rubber  (0.068-in.), 
Teflon  (0.033-in),  and  polyethylene 
plastic  (0.007-in.)  - were  evaluated  as 


friction-reducers  in  a series  of  tests  on 
saline-ice  specimens  at  -IOC.  A constant 
strain  rate  of  8.5  x 10'^  in./in.  sec'^  was 
maintained  throughout  the  tests.  The 
results  are  also  shown  in  Figures  14  and 
15.  While  these  tests  produced  the 
expected  result  of  lowering  the  compres- 
sive strength,  shear  planes  were  still  the 
dominating  failure  mechanism  of  hori- 
zontal ice  specimens  with  vertical 
cleavage  cracking  notably  absent.  How- 
ever, lateral  bulging  of  most  horizontal 
specimens  was  reduced,  and  this  would 
probably  have  a small  effect  on  longi- 
tudinal deflection  measured  by  an 
extensometer  attached  to  the  surface  of 
the  ice  sample.  It  was  observed  that  the 
thicker  rubber  caused  plastic  flow  in 
several  ice  specimens  even  at  this 
moderate  strain  rate,  while  cardboard 
absorbed  moisture  and  had  to  be  replaced 
often.  It  is  recommended  that  a thin 
sheet  of  some  friction-reducing  material 
(such  as  plastic  or  Teflon)  be  used  during 
compression  testing  of  sea  ice;  however, 
previously  mentioned  compressive 
strength  tests  without  friction-reducers 
are  still  valid,  since  they  provide  an 
upper  limit  for  design  considerations  of 
ice  action  on  a structure. 


Flexural  Strength 

Both  laboratory  and  field  tests  were 
performed  to  determine  the  flexural 
strength  and  apparent  elastic  modulus  of 
sea-ice  beams  (1,  3,  4,  12).  Laboratory 
specimens  were  all  2 x 2 x 16-in.  span 
beams,  tested  under  a symmetrical  two- 
point  load,  spaced  6 in.  apart.  The  test 
setup  is  shown  in  Figure  18.  All  beams 
were  sawed  from  the  ice  block  with  their 
longitudinal  axis  parallel  to  the  ice 
surface  (horizontal  beams),  and  tested  in 
an  orientation  that  placed  the  normal 
vertical  ice  crystal  growth  in  the  same 
plane  as  the  applied  load.  Once  roughly 
cut  from  the  ice  block,  these  beams  were 


finished  on  a planer  to  a nominal  2x2- 
in.  cross  section.  However,  analysis  was 
based  on  calculating  the  section  modulus 
of  each  beam  from  its  actual  cross- 
sectional  dimensions.  Five  test  tempera- 
tures (-4,  -7,  -10,  -20,  and  -27C)  and  one 
salinity  level  (7  to  9 ppt)  were  the 
independent  variables  for  all  laboratory 
flexural  strength  tests.  A strain-gage- 
type  extensometer  with  a 4.5-in.  gage 
length  attached  to  the  underside  of  the 
beam  provided  continuous  strain  output. 
The  gage  was  positioned  symmetrically 
with  the  beam  centerline  to  account  for 
only  strain  within  the  constant  moment 
section  of  the  beam. 

Flexural  strength  data  on  annual  sea 
ice  at  McMurdo,  Antarctica,  and  Barrow, 
Alaska,  come  from  testing  large  in-situ 
beams  approximately  40  in.  wide  by  66  to 
96  in.  deep  and  ranging  from  50  to  100  ft 
in  length.  The  apparatus  for  cutting  and 
testing  large  in-situ  beams  is  shown  in 
Figures  19  and  20.  The  ice  saw  (Figure 
19)  is  capable  of  cutting  ice  up  to  10  ft 
thick,  and  the  load  frame  (Figure  20)  can 
apply  loads  up  to  200,000  lb,  which  is 
measured  with  a strain  gage  load  cell.  To 
avoid  the  necessity  of  including  shear 
deformation  in  the  analysis,  all  of  the 
large  field  beams  maintained  a length- 
to-depth  ratio  between  8 and  10.  These 
beams  were  tested  under  three  boundary 
conditions  (cantilever,  fixed-fixed,  and 
simply  supported)  with  an  LVDT  having  a 
6-in.  stroke  that  continually  monitored 
deflections.  Except  for  a few  beams 
tested  at  -2.5C  (essentially  isothermal), 
the  remainder  of  the  large  field  beam 
tests  were  conducted  with  a thermal 
gradient  across  the  beam.  Ice-sheet 
temperature  as  well  as  salinity  profiles 
were  continuously  logged  during  the 
tests.  The  general  salinity  data  range  for 
all  field  tests  was  5 to  7 ppt. 

In  addition  to  the  large  in-situ 
beams,  several  medium-sized  beams  with 
a cross  section  from  30  to  60  sq  in.  were 
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T«t  apparatus  for  dttarmining  fitxural  itrenflth  and  apparent  elastic 
modulus  of  2 X 2 X IS'in.  span  beam. 


extracted  from  the  ice  sheet  and  center- 
point  loaded  as  simply  supported  beams. 
Several  groups  of  small  2 x 2 x 16-in. 
TCams  were  also  cut  from  ice  sheet 
bl^Ks  and  tested  identically  to 
laboratory  experiments.  The  typical 
stress  rate  for  applying  the  failure 
flexural  load  in  the  field  was  greater 
than  1<»  psi/sec,  while  typical  stress  rates 
for  laboratory  beam  tests  ranged  from  25 

* iM  Sea  iuc  lu  5u 

psi/sec  for  cold  (-27C)  sea  ice. 

The  flexural  strength  from  both 
laboratory  and  field  tests  is  shown  in 
Figure  21  as  a function  of  the  square 
root  of  the  brine  volume  present  in  the 
ice.  The  same  analytical  technique  used 
to  relate  tensile  and  compressive 
strength  to  tf»e  combined  effect  of 
salinity  and  temperature  was  employed. 


The  data  were  plotted  using  the 
individual  field  data  points  and  the  mean 
value  for  each  laboratory  sample  group. 
A weighting  factor  was  applied  to  give 
more  significance  to  the  field  tests  by 
assuming  25  laboratory  specimens 
equaled  one  large-scale  field  test.  The 
eqi;ation  of  the  straight  line  based  on 
least-square  analysis  is: 

Oj  = 139.1-8.82  7^  (5) 

where  v is  the  brine  volume  in  ppt  and 
oj  is  the  flexural  strength  in  Ib/sq  in. 
Again,  each  isothermal  laboratory  data 
poup  is  shown  with  a 95%  confidence 
level  about  its  mean.  If  one  extrapolates 
the  best  fit  line,  it  is  found  that  Equation 
5 predicts  a flexural  strength  (modulus  of 
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Figof#  20.  Hydraulic  test  fr«na  with  200.(X»lb  capadty  for  loading  large  in  titu  lerice  beam. 
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rupture)  of  about  140  psi  for  ice  of  zero 
brine  content  but  the  same 
microstructure  as  sea  ice. 

Figure  22  shows  the  effect  of 
temperature-salinity  on  both  flexural  and 
horizontal  tensile  strength  as  a function 
of  the  square  root  of  brine  volume.  It  is 
significant  that  the  analysis  for  each 
case  resulted  in  first  order  equations 
that  have  nearly  the  same  line  slope, 
that  is,  -8.82  and  -9.99.  Of  course,  this  is 
as  it  should  be  since  the  portion  of  the 
beam  below  the  neutral  axis  is  under 
tension  and  represents  a tension  failure. 
The  indicated  ultimate  tensile  strength 
of  the  ice  based  on  beam  flexural 
formula  is  approximately  20  psi  greater 
than  the  observed  uniaxial  tensile 
strength  (horizontal).  This  disparity  is  a 
common  characteristic  of  brittle 
materials  that  do  not  have  the  same 
tension  and  compression  properties.  The 
neutral  axis  in  the  beam  specimen  must 
shift  towards  the  stiffer  side  of  the  beam 
(compression  side  for  ice)  in  order  to 
maintain  equality  of  the  resultants  of  the 
tension  and  compression  forces.  Thus, 
the  fiber  stress  computed  by  the  flexural 
formula  is  less  than  the  true  fiber  stress 
on  the  stiffer  side  of  the  beam  and  is 
greater  than  the  actual  fiber  stress  on 
the  other  side. 

Elastic  Modulus 

In  addition  to  strength  properties, 
the  single  most  important  material 
property  of  sea  ice  is  its  clastic  (Young's) 
modulus,  which  is  necessary  for 
analyzing  ice  plate  behavior.  Conse- 
quently, strain  values  obtained  during  the 
beam  tests  were  used  to  determine  the 
apparent  modulus  of  elasticity  for  sea 
ice  as  a function  of  brine  volume.  All 
laboratory  strain  data  were  measured 
with  extensometers  bonded  tc  the  under- 
side of  2 X 2 X 16-in.  simply  supported 
beams  subjected  to  two-p>oint  loadings.  In 


addition  to  the  laboratory  experiments, 
the  large-scale  field  beam  tests 
measured  deflections  with  a 6-in.  stroke 
LVDT.  These  results  are  displayed  in 
Figure  23,  and  a straight-line  approxima- 
tion using  least  squares  is  calculated  for 
all  flexural  tests. 


771  - 63.2 


where  Ef  is  the  apparent  flexural  modu- 
lus of  elasticity  in  Ib/sq  in.  and  v is  th#* 
brine  volume  in  ppt. 

Another  method  for  finding  the 
modulus  of  elasticity  is  to  measure  the 
strain  corresponding  to  a given  streis 
state  in  a right-circular  cylinder  subjec- 
ted to  either  uniaxial  compression  or 
tension.  During  compressive  strength 
testing,  deflection  measurements  were 
recorded  with  extensometers  having  a 
gage  length  of  2.4  in.  that  were  attached 
directly  to  each  specimen  by  bonding 
with  freshwater  ice.  Hov/ever,  elastic 
modulus  values  found  by  compression 
testing  are  strain-rate  sensitive.  The 
results  of  three  different  compressive 
stiain-rate  tests  performed  at  one  brine 
volume  value  on  horizontal  sea-ice  speci- 
mens are  shown  in  Figure  23.  The  moduli 


averages  for  a limited  statistical  sample 
for  both  1.7  x 10*’  and  3.4  x 10*^  in. /in. 
sec**  strain  rates  are  shown  as  triangles 
on  Figure  23.  Since  these  values  bracket 
the  flexural  data,  these  two  strain  rates 
are  deserving  of  a more  complete  evalu- 
ation. Another  possible  contribution  to 
greater  compressive  deformation 
measurement  in  horizontal  specimens 
that  cause  unrealistically  low  elastic 
moduli  values  is  end  constraints.  No 
investig=»tions  were  made  of  end-con- 
straint effects  on  the  elastic  modulus  of 
sea  ice.  in  addition  to  further  strain-rate 
and  end-constraint  testing,  better  strain 
and  deflection  measurement  techniques, 
especially  at  warm  temperatures,  are 
required  to  properly  assess  the  deforma- 
tion behavior  of  sea  ice. 
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Figure  21 . Flexural  strength  versus  brine  volume. 
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Figure  22.  Flexural  ttreitgth  (rupture  modulus)  and  horizontal  tensile  strength 
versus  brine  volume. 
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Figure  23.  Young's  modulus  of  sea  ice  versus  brine  volume. 


Figure  24.  Laboratory  setup  for  plate  experiments. 


Several  "^ea-ice  plate  experiments 
were  performed  in  the  laboratory  to 
investigate  the  behavior  of  a CEL- 
developed  sea-ice  strain  transducer  (lf>) 
and  to  determine  the  modulus  of  elastic- 
ity by  monitoring  strains  and  deflections 
of  sea-ice  plates  under  load.  Simply 
supported  sea-ice  plates  (9  ppt  salinity), 
4. .*5  in.  thick  and  42  in.  in  diameter,  were 
subjected  to  a circular  ring  loading 
applied  by  a hydraulic  ram.  The 
laboratory  setup  is  shown  in  Figure  24, 
while  Figure  25  illustrates  the  geometry, 
ring  loading,  and  location  of  embedded 
strain  gages  common  to  all  plate  tests. 
The  gages  were  grouped  in  two  pairs  so 
that  each  pair  measured  strain  in  two 
orthogonal  directions  to  cancel  any 
anomalies  arising  from  asymmetry. 
Gages  1 and  4 were  made  of  two  active 
arms  that  were  temperature- 
compensated  and  encased  in  a pliant 


material  called  Gagekote  5,  while  gages 
2 and  3 were  identical  except  that  the/ 
were  encased  in  Plexiglas.  All  four  gages 
had  an  active  length  of  5.5  in.  and  a gage 
factor  of  2.0.  Their  development  and 
testing  are  discussed  in  greater  detail  in 
Reference  14. 

In  order  to  induce  elastic  behavior, 
loading  was  applied  at  stress  rates 
ranging  from  7.5  to  15  psi/sec  until  a 
preselected  peak  load  was  reached;  then 
the  load  was  released.  The  strains  were 
continuously  recorded,  and  peak  strains 
were  easily  distinguished.  This  pattern 
was  repeated  for  different  peak  loads. 
Some  scatter  of  peak  load  versus  peak 
strain  was  noted  for  all  gages  during 
different  tests.  However,  the  top  gages 
in  the  compression  zone  exhibited  a 
linear  relationship  with  the  load,  whereas 
the  tension  gages  suggested  a nonlinear 
response,  perhaps  due  to  tension  cracking 
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of  the  sea  ice.  In  addition  to  the  embed- 
ded strain  gages,  deflections  were 
monitored  by  a single  centerline  LVDT 
with  a 6-in.  stroke  plus  two  extenso- 
meters  with  gage  lengths  of  2.5  in. 
bonded  to  both  the  top  and  bottom  of  the 
sea-ice  plates. 

Plate  tests  were  performed  at  four 
constant  temperatures:  -4,  -10,  -20,  and 
-30C.  The  elastic  modulus  values  for 
each  test  temperature,  which  were  cal- 
culated from  plate  theory  equations,  are 
shown  in  Figure  26.  The  averaged  values 
from  gages  1 and  k provided  the  strain 
input,  while  the  output  from  the  LVDT 
and  both  extensometers  were  averaged 
to  supply  deflection  data.  The  apparent 
elastic  moduli  found  from  the  flexural 
tests  are  superimposed  on  this  same 
figure  for  comparison.  In  a later  section 
on  analytical  methods,  output  from  the 
elastic  finite-element  computer  code 
will  be  compared  with  the  sea-ice  plate 
strain  data. 

FractiM’e  Toughness 


either  the  stress  intensity  factor,  K,  or 
the  strain  energy  release  rate,  y , are 
known,  the  stress  field  at  the  crack  tip 
becomes  completely  defined.  The  factor 
K increases  with  an  increase  in  applied 
load,  but  a small  area  of  plastic  de- 
formation ahead  of  any  growing  crack 
can  partially  offset  the  intensification  of 
the  local  stress  field  caused  by  the 
applied  load,  P,  and  crack  length,  c. 
However,  K increases  faster  than  the 
plastic  strain  area,  and  its  value  rises 
rapidly  during  this  initial  period  of  slow 
crack  growth.  At  some  combination  of  P 
and  c,  represented  by  a critical  value  of 
K,  the  stress  intensity  at  the  crack  edge 
becomes  so  large  that  the  crack  can 
propagate  without  supplying  additional 
energy.  At  this  time  P reaches  a maxi- 
mum, the  crack  becomes  unstable  and 
propagates  rapidly  instead  of  slowly.  The 
critical  value  of  K is  called  the  fracture 
toughness,  Kj  , of  the  material  and  is 
related  to  the  critical  strain  energy 
by 


In  the  analysis  of  ice  sheet  behavior 
under  sudden  or  sustained  loading,  an 
allowable  flexural  stress  has  been  tradi- 
tionally considered  as  the  design  (or 
application)  criterion.  Normally,  this 
allowable  stress  condition  represents  an 
arbitrary  reduction  in  a known  ultimate 
flexural  strength,  usually  an  average 
value  of  a set  of  laboratory  and  field 
experiments.  Both  fracture  mechanisms 
and  failure  modes  of  sea  ice  are  being 
studied  in  an  attempt  to  define  a suitable 
failure  criterion  for  ice  sheets.  Once  a 
failure  criterion  is  determined,  a factor 
of  safety  can  be  established  with  greater 
confidence. 

One  possible  approach  to  the  ice- 
sheet  cracking  problem  is  to  treat  it  by 
Griffith  theory  as  a brittle  fracture  (13). 
First,  the  point  is  defined*  at  which  an 
existing  crack  becomes  unstable.  Once 


(7) 


where  E is  Young's  modulus  and  p is 
Poisson's  ratio.  According  to  Griffith, 
cracks  will  propagate  rapidly  and  pro- 
duce a brittle  fracture  (the  area  of 
plastic  strain  is  very  small)  when  the 
applied  stress  reaches  a failure  value  of 
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or  rewritten  using  Equation  7 
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(8b) 
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t Figure  26.  Sea  ice  elastic  modulus  versus  ice  temperature. 


where  c is  the  critical  crack  length.  surface.  A binomial  expression  can  be 

While  Yc  difficult  to  evaluate,  written  to  calculate  the  fracture 

I K|j,  , the  fracture  toughness  property,  is  toughness  as  a function  of  the  notch-to- 

easier  to  find  by  experimentation.  To  beam-depth  ratio  (15). 
determine  K|  for  sea  ice  as  a function 
of  brine  volume,  laboratory  tests  were 

performed  on  2 x 2 x 16-in.  beams  with  3 P b ' ^ ( £..  - 1^) 

central  saw-cut  notches  0.5  in.  deep.  A K,  = =-! —6  (9) 

I two-point  loading  with  a 6-in.  spacing  c 2 w d^ 

\ was  applied  at  a loadhead  speed  of  0.2 

in./min.  Approxi'nately  25  notched  sea-  where  P is  the  applied  load,  b is  the 
ice  beams  were  tested  at  each  of  two  notch  depth,  I j is  the  beam  suoport 
test  temperatures:  -20  and  -IOC.  All  of  spacing,  1 2 is  the  two-point  load 
the  beams  failed  instantly  with  a crack  spacing,  w is  the  beam  width,  d is  the 
initiating  at  the  root  of  the  notch  and  beam  depth,  and 
rapidly  propagating  through  to  the  top 
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Figure  27.  Fracture  toughness  versus  brine  volume. 


B = 1.99  - 2.<»7(b/d)  + 12.97(b/d)^ 

- 23.17(b/d)^  + 2<».80(b/d)'' 

= 1.918  (for  b/d  = 0.25)  (10) 

Using  the  geometric  constants  from  the 
insert  of  Figure  27  and  Equation  10, 
Equation  9 can  be  rewritten  as 

K,  = 2.5<»3P  (11) 

c 

The  fracture  tqughness  results  from 
the  laboratory  experimentation  are 
shown  in  37  as  a function  of  brine 

volume,  V , The  straight-line  approxima- 
tion is  calculated  by  the  linear  regression 


technique  of  least  squares  and  represents 
a tentative  relationship  between 
and  V. 

Kj  = 129-12.<»  (12) 

c 

w^here  K is  the  material  property 
fracture  toughness  expressed  in  psi  >/Tn. 
and  V is  the  brine  volume  in  ppt.  Using 
Figure  26  to  obtain  a range  for  the 
modulus  of  elasticity  and  letting  u = 0.3, 
one  can  calculate  the  critical  energy  Yc 
which  causes  a crack  to  propagate  with- 
out any  additional  energy  required,  from 
Equation  7.  For  sea  ice,  » 1 x 10*2 
in.-lb/sq  in.  The  size  of  a critical  crack 
can  be  found  by  rewriting  Equation  8b  as 
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c = V (13) 

TT  O 

If  nominal  values  of  65  psi  VitiT  and  100 
psi  are  chosen  for  K|^  and  o respec- 
tively, then  the  critical  crack  length 
would  be  approximately  0.25  in.  Of 
course,  sea  ice  contains  brine  inclusions 
that  are  around  0.05  in.  in  diameter, 
considerably  smaller  than  the  0.25  in. 
There  is  no  doubt  that  these  brine 
pockets  serve  as  stress  concentration 
points  and  provide  a source  for  crack 
initiation.  However,  it  does  seem 
probable  that  some  dislocation  move- 
ment due  to  localized  plastic  strain  could 
occur  that  would  cause  slow  crack  propa- 
gation until  the  critical  crack  length  is 
reached  by  linking  several  brine  pockets 
with  small  subcritical  cracks. 

It  is  not  sufficient  to  attempt  to 
formulate  a failure  criterion  for  sea-ice 
sheets  solely  on  the  basis  of  fracture 
mechanics  theory  and  small  notched- 
beam  tests  in  the  laboratory.  As  in  other 
material  properties  experiments,  how- 
ever, it  is  necessary  to  investigate  and 
understand  microscopic  behavior  before 
it  is  applied  to  the  macroscopic  scale  of 
a floating  ice  sheet. 

Creep  Behavior 

It  is  well  known  that  the  assumption 
of  clastic  behavior  for  sea  ice  is  valid 
only  for  a narrow  range  of  load  magni- 
tude and  length  of  application.  Many 
previous  studies  have  considered  visco- 
elastic ice  behavior  (16,  17,  18).  At  CEL 
two  laboratory  investigations  have  been 
made  on  sea-ice  creep  response:  one  on 
simply  supported  beams  (3),  and  one  on 
uniaxial  compressive  and  tensile 
specimens  (19).  In  addition,  field  exper- 
imentation on  long-term  sea-ice  behavior 


of  in-situ  cantilever  beams  has  been 
performed  as  a large-scale  verification 
of  laboratory  efforts. 

Uniaxial  Laboratory  Tests.  An 
extensive  laboratory  experimental  pro- 
gram has  been  performed  to  define 
uniaxial  deformation  time  history  of  sea 
ice.  Both  uniaxial  compressive  and 
tensile  creep  tests  were  conducted  on 
cylindrical  specimens  having  a length-to- 
depth  ratio  of  2:1  that  were  cored  either 
horizontally  or  vertically  from  sea-ice 
blocks.  A constant  load  was  applied  to 
each  specimen  through  a deadweight 
lever  arm  system  (commonly  used  in  soil 
compaction  testing).  This  applied 
constant  load  was  predetermined  as  a 
percentage  of  the  average  elastic 
compressive  or  tensile  strength  for  a 
given  ice  temperature  and  crystal  orien- 
tation. Three  different  percentages  (10, 
25,  50%)  were  considered  as  load  level 
parameters.  Several  tests  at  each 
parametric  combination  of  loading  condi- 
tions, crystal  orientation,  and  load  levels 
were  conducted  as  a series  for  each  of 
four  ice  temperatures  (-4,-10,-20,-270. 
Consequently,  a total  of  48  differ- 
ent parametric  combinations  had  to  be 
considered. 

Ice  growth  has  been  controlled  to 
maintain  a constant  salinity  range, 
density,  and  crystal  size  representative 
of  natural  sea  ice.  Still,  for  time-depen- 
dent tests  this  became  an  awesome  task 
requiring  multiple-specimen  testing  over 
a single  time  interval.  Two  test  frame 
setups  (Figures  28  and  29)  were  devised 
to  accommodate  four  compression  and 
two  tension  tests  at  the  same  time.  Doth 
frames  are  deadweight  systems  with  the 
compression  loading  in  Figure  28  applied 
through  a 40:1  lever  arm.  Figure  29 
shows  a 10:1  lever-arm  compression 
loading  frame.  Still,  over  300  specimens 
had  to  be  tested  to  assure  piroper 
statistical  sampling. 
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Figure  30.  Viscoelastic  behavior  of  sea-ice  specimens  with  compression  loading 
and  horizontal  (perpendicular)  crystal. 


The  output  for  each  individual  test  10  and  25%  stress  levels,  since  the  50% 
is  the  strain  history,  shown  in  Figures  30  stress  level  produced  failure  within  5-10 
through  33,  but  each  curve  represents  an  min.  On  the  other  hand,  the  tensile  strain 
average  for  a sample  size  ranging  from  3 histories  in  Figures  32  and  33  are  shown 
to  6 specimens.  The  deflection  is  mea-  for  those  stress  levels  that  produced 
sured  by  extensometers  with  a gage  good  creep  data.  First,  a 10%  stress  level 
length  of  2.k  in.,  attached  directly  to  the  produced  little  creep  response  in  tension 
ice  specimens.  Each  figure  displays  specimens  over  the  30-hr  test  period,  so 
different  temperature  and  stress  level  several  attempts  were  made  to  substi- 
curves  for  a specified  loading  condition  tute  a 75%  stress-level  loading.  How- 
and  crystal  orientation.  ever,  most  75%  stress-level  loadings 

Theoretically,  each  figure  should  produced  instantaneous  failure.  Several 
display  12  curves,  but  higher  stress  levels  50%  stress-level  curves  are  also  missing 
fail  ice  specimens  very  quickly.  There-  due  to  instantaneous  failure,  while 
fore,  those  fast-acting  and  virtually  absent  25%  curves,  especially  at  lower 
nonexistent  strain  histories  have  been  temperatures,  represent  no  observable 
deleted.  In  Figures  30  and  31  compres-  response  over  the  testing  time  period.  It 
sive  strain  histories  are  shown  for  only  is  intended  to  fill  out  data-deficient 
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Figure  31.  Viscoelastic  behavior  of  sea-ice  specimens  with  compression  loading 
and  vertical  (parallel)  crystal. 


tension  tests  if  required  to  adequately 
predict  time-dependent  ice-sheet  per- 
formance. 

While  conducting  uniaxial  compres- 
sion tests  on  sea  ice,  one  can  observe  the 
failure  planes  develop  as  a function  of 
time.  Generally,  specimens  loaded 
parallel  to  their  preferred  crystal  orien- 
tation failed  along  vertical  cleavage 
planes,  but  macroscopic  fracture 
occurred  only  at  25  and  50%  stress 
levels,  while  the  10%  stress  level 
produced  only  microscopic  cracking.  The 
horizontal  specimens,  loaded  perpendic- 
ular to  ice  crystal  growth,  universally 
failed  along  developed  shear  planes 
independent  of  ice  temperature  and 


stress  level.  However,  different  combi- 
nations of  these  two  parameters  did  have 
a significant  effect  on  the  rapidity  of 
failure. 

Laboratory  Beam  Tests.  Observa- 
tions were  concerned  with  both  creep 
deflection  and  strain  measurements  of  2 
X 2 X 18-in.  sea-ice  beams  with  a two- 
point  load  applied  at  the  one-third 
point* , However,  only  the  strain  mea- 
surements over  a 2.4-in.  gage  length  will 
be  considered  in  this  discussion.  Strain- 
gage-type  extensometers  were  used  for 
continuous  strain  monitoring  during 
testing. 

To  determine  creep  properties  it 
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Figure  32.  ViKoelastic  behavior  of  tea-ice  specimens  with  tension  loading  and 
horizontal  (perpendicular)  crystal. 


was  found  n cessary  to  enclose  the  test 
apparatus  holding  the  beam  in  a plastic 
shroud  (1)  to  shield  the  beam  and  the 
instrumentation  from  the  effects  of  air 
movement  within  the  test  chamber  and 
cyclic  temperature  variation  {*  2.2C)  of 
the  chamber,  and  (2)  to  minimize 
sublimation  of  the  specimen  during  the 
long  term  of  the  test.  Continuous  tem- 
perature monitoring  within  the  plastic 
enclosure  indicated  fluctuation  was 
reduced  to  a nominal  ♦1/2C.  The  test 
apparatus  for  the  1 8-in.-span  beam  was  a 
small  soil-test  frame  shown  in  Figure  34 
which  had  been  modified  to  accommo- 
date a beam  specimen.  The  end  supports 
for  the  beam  and  the  loading  yoke, 
having  multiaxis  freedom  of  movement 


for  alignment,  had  rounded  bearing 
surfaces  to  minimize  the  resistance  of 
the  specimen  to  rotation.  Thin  metal 
shims  were  inserted  between  the  ice 
beam  and  load  bearing  points  to  prevent 
pressure  penetration  of  the  ice. 

Creep  strain  histories  are  presented 
in  Figure  35  for  two  stress  and  tempera- 
ture parameters.  Each  curve  represents 
an  average  of  five  individual  beam  tests. 
The  salinity  range  was  7 to  9 ppt  for  all 
beam  creep  tests  measuring  strain. 

Cracks  in  the  tension  face  of  the 
beam  occurri.ng  particularly  at  crystal 
boundaries  were  frequently  found  during 
inspection  of  the  specimen  upon  test 
completion.  An  effective  method  of 
detecting  the  surface  and  internal  cracks 


10 


o • 40  p$i  150%)  »t  -4C 


Figur«  33.  ViKoelastic  behavior  of  wa-ice  ipecimcns  with  ttntion  loading  and 
vertical  (parallell  crystal. 


at  the  time  of  occurrence  has  yet  to  be 
developed;  however,  even  with  such  a 
capability,  analysis  of  the  stress<strain 
state  resulting  from  crack  development 
would  be  difficult.  A rarely  observed 
creep  phenomenon,  that  of  deformation 
and  slippage  in  the  crystal  lattices 
occurring  at  a triple  crystal  intersection, 
is  shown  in  Figure  36. 

Field  Beam  Tests.  This  large-scale 
experiment  was  designed  to  provide  data 
jXJints  to  help  verify  the  extensive 
laboratory  program  defining  the  visco- 
elastic properties  of  sea  ice.  Tests  were 
conducted  on  large  in-situ  cantilever  ice 
beams  cut  from  annual  sea  ice  - 


ungrounded,  level,  and  generally  crack- 
free.  Beams  had  cut  dimensions  of 
approximately  3 ft  in  width  and  ft 
long.  The  ice-sheet  thickness  was  75  in. 
during  the  testing  period.  A 30-hp, 
ladder-type  trencher  with  specially 
designed  conical  teeth  cut  the  beams, 
leaving  8-in.-wide  open  water  trenches. 
However,  it  became  necessary  to  employ 
heat-tape  floating  in  the  trenches  to 
prevent  freezeback. 

A deadweight,  <»:1  lever  arm  system 
provided  an  instantaneously  applied, 
constant,  downward  load  to  the  free  end 
of  the  beams.  The  actual  load  applied  to 
each  beam  was  continuously  recorded  by 
a load  cell,  while  beam  deflection  was 
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Figure  34.  Small-beim  creep  test  freme  ttrein-gage-type  extensometer 
measures  strain  over  a 2’in.  gage  length. 

monitored  at  the  load  application  point  (Figure  38  and  39,  respectively)  were 
by  an  LVDT  (linearly  varying  displace*  displaying  elastic  behavior  after  the 
ment  transducer).  The  entire  setup  is  second  step  loads  at  stress  levels  of  27.9 
shown  in  Figure  37.  The  data  are  plotted  and  30.6  psi  applied  at  times  45  and  65  hr 
in  Figures  38  and  39.  A salinity  profile  into  the  tests,  respectively.  On  the  other 


was  recorded  for  each  of  the  seven  ice  hand,  beams  2 and  6 (Figure  38)  showed 
beams  with  salinities  falling  in  a r^arrow  slightly  plastic  behavior  under  their  third 
range  between  5.2  to  5.7  ppt.  Ice  surface  step  loads  occurring  at  23  and  40  hr  and 


temperatures  ranged  from  -12  to  -20C  at  stress  levels  of  23.9  and  42.8  psi, 
during  the  test  period.  During  each  test  a respectively.  This  very  limited  informa- 
temperature  profile  of  the  ice  sheet  was  tion  is  not  sufficient  to  adequately 
periodically  recorded.  define  elastic-plastic  beltavior  twunda- 

Each  test  was  run  for  a period  of  45  ries.  However,  it  appears  that  for  sea  ice 
to  70  hr,  with  step  loads  added  at  at  temperatures  between  -10  and  -I5C 


I different  times  in  an  attempt  to  identify  inelastic  behavior  can  occur  at  relatively 

[ boundaries  of  viscoelastic-viscoplastic  low  stress  levels  applied  over  a time 

I behavior  as  a function  of  stress-  period  greater  than  24  hr.  A scries  of 

temperature  levels.  All  the  beams  controlled  laboratory  tests  should  be 

showed  an  elastic  response  after  their  performed  to  determine  viscoelastic  and 

initial  step  loads,  and  both  beams  4 and  7 plastic  limits  on  sca-ice  behavior. 
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Figure  37.  Test  setup  for  field  effort  to  define  viscoelastic  sea  ice  behavior. 
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Figure  38.  Creep  time  history  for  in-sito  cantilever  beams  2. 4.  and  6. 
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Figure  39.  Creep  time  history  for  in-situ  cantilever  beams  3,  5,  and  7. 


Of  course,  field  tests  cannot  have 
the  controlled  temperature  environment 
of  the  laboratory;  in  fact,  when  a beam 
is  cut  in  situ,  the  sides  are  exposed  to 
seawater,  slowly  warming  the  beam  from 
three  directions.  Consequently,  the  ice 
beam  is  continually  warmed  throughout 
the  test;  therefore,  in  order  to  discuss 
stress-temperature  relationships,  an 
average  temperature  gradient  was 
determined  for  each  test. 

The  creep  curves  (deflect ion- time 
history  of  the  six  cantilever  test  beams) 
can  be  given  in  the  form  of  strain 
functions  based  solely  on  beam  geome- 
tries. First,  a displacement  function, 
u(x),  is  assumed  to  be  of  the  form 


2 

u(x)  = A + A.x  + A^x 
o I 2 

+ AjX^  (14) 


and  the  coefficients  can  be  found  from 
evaluating  boundary  conditions.  Then, 
strain  can  be  written  as 

c(,,,)  y (15) 

3x^ 
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Figure  40.  Loading  geometry  for  field  testing  of  viscoelastic  sea-ice  beams. 


where  y is  the  distance  from  the  extreme 
beam  fiber  to  its  neutral  axis.  For  the 
field  ice  beams  with  the  geometric  con- 
figuration of  Figure  ^0,  the  strain  func- 
tion becomes 

cl,,,)  , _Av  <x  - a)  A(tl_  ,1^, 

(H-cr  (21*  c-  3a) 


an  analytical  approach  based  on  two- 
dimensional  elasticity  which  can  easily 
handle  the  temperature  gradient  across 
the  ice-sheet  thickness.  More  recently, 
the  finite-element  technique  has  been 
expanded  to  include  viscoelastic,  as  well 
as  elastic,  behavior  of  sea-ice  sheets. 
All  three  analytical  approaches  arc 
summarized  in  the  following  sections. 


where  x > a,  H is  the  beam  length,  a and 
c are  the  distances  from  the  free  end  to 
the  points  of  load  application  and  deflec- 
tion measurement,  respectively,  and  A(t) 
is  the  deflection  time  history  for  each 
individual  beam  test.  Equation  16  will  be 
used  in  a later  section  on  viscoelastic 
finite-element  analysis  to  determine 
relaxation  functions  required  for  input  to 
the  computer  code. 


ANALYSIS  OF  FLOATING 
ICE  SHEETS 

Initially,  the  problem  of  safe  bearing 
capacity  of  floating  ice  sheets  was 
analyzed  by  thin-plate  theory,  assuming 
the  it  e to  be  an  elastic  plate  resting  on 
an  elastic  foundation  of  the  Winkler 
type.  Then,  development  of  an  elastic 
finite-element  computer  code  provided 


Elastic  Plate  on  an  Elastic 
FourKlation 

A floating  annual  ice  sheet  of 
infinite  size  can  be  analyzed  by  elastic 
theory  as  a plate  resting  on  an  elastic 
foundation  (20,  21).  The  ice  sheet  has 
flexural  rigidity;  as  it  deflects,  there  is  a 
restoring  force  due  to  the  water  pressure 
that  is  proportional  to  the  deflection. 
Simplifying  assumptions  in  the  criteria 
were;  (I)  the  ice  behaves  as  an  isotropic 
medium,  (2)  its  relaxation  characteristics 
are  outside  the  domain  of  the  loading 
times,  and  (3)  the  moving  surface  loads 
do  not  create  resonance  wave  conditions 
in  the  plate. 

Plate  theory  with  small  deflections 
may  be  used  to  determine  the  critical 
stress  causing  the  initial  cracking  of  the 
plate.  Hertz  (22)  was  the  first  to  solve 
this  problem  for  a uniform  load 
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distributed  over  a circular  area.  Wyman 
(23)  expressed  the  solution  in  terms  of 


determined  numerically  that  when  a/h  > 
1.724,  the  critical  stress  given  by  thin- 


Kelvin  functions,  where  +he  maximum  plate  theory  is  valid.  For  a/h  < 1.724,  the 


stress  and  deflection  occurring  under  the 
loading  center  are 


1 + (a/6)  ker'  (a/6) 


radius  of  the  load  circle  should  be 
recomputed  by  the  following  equation. 


aj  = (1.6a^  + h^) 


0.675h 


3P(1  +n)  kei-  (a/6) 
„k2  (a/6) 


where  aj  = adjusted  load  circle  radius 

a = actual  load  circle  radius 
h = ice  sheet  thickness 


where  w = the  maximum  vertical 

max  , 

deflection 

o„  = the  maximum  tensile 
max 

stress 

P = the  total  load 

k = the  foundation  modulus 

(unit  weight  of  seawater) 

a = the  radius  of  load 

distribution 


[12k  (1  J 

characteristic  length 
called  action  radius 

= the  plate  thickness 

= the  modulus  of  elasticity 

- Poisson's  ratio 


Of  course,  superimposing  the  load 
effects  of  several  aircraft  landing  gear, 
vehicle  wheels,  or  tracks  is  necessary  to 
determine  the  actual  stress  and  deflec- 
tion in  the  ice  sheet.  The  superposition 
of  the  components  of  stress  for  the 
remaining  loads  onto  the  initial  stress 
condition  produced  by  the  load  at  the 
zero  coordinate  point  and  given  by  Equa- 
tion 18  can  be  accomplished  by  a summa- 
tion of  partial  derivatives  of  the  bending 
effect  relative  to  the  coordinate  position 
of  the  particular  load.  If  we  let  vj  = 
Wj/P,  where  Wj  = plate  deflection 
produced  by  any  one  of  the  off-zero 
coordinate  loads,  then 


6x^  ~ 6x^  6y^  ^ 6y^  ^ 


One  limitation  on  the  above  theory 
is  that  the  radius,  a,  of  the  loaded  area 
should  be  large  in  comparison  to  the 
thickness  h of  the  plate.  Tracked 
vehicles  usually  meet  this  requirement, 
but  wheeled  vehicles  or  aircraft  do  not. 
For  a small  radius  of  loading,  thin-plate 
theory  predicts  a stress  that  is  too  high. 
Westergaard  (24)  investigated  this 
problem  for  a simply  supported  plate  and 


The  total  deflection  of  the  ice  plate 
from  loading  by  aircraft  or  vehicles  was 
calculated  by  the  same  general  technique 
used  for  stress  calculation,  starting  with 
Equation  17  for  the  deflection  of  the 
zero-coordinate  load. 
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Elastic  Finite  Element  Method 

To  avoid  the  limiting  assumption  of 
thin-plate  theory  with  regard  to  the 
loading  circle  radius,  to  incorporate 
thermal  gradient  effects  of  the  ice 
sheet,  and  to  provide  a versatile 
numerical  tool  for  evaluating  ice-sheet 
bearing  capacities,  an  elastic  finite- 
element  computer  code  was  developed  at 
CEL.  The  program  is  based  on  two- 
dimensional  elasticity  with  the  structural 
system  idealized  as  an  axisymmetric 
solid  composed  of  an  arbitrary  number  of 
layers  (4,  25).  These  layers  may  have 
orthotropic  or  temperature-dependent 
material  properties.  The  solid  is  sup- 
ported by  a fluid  foundation  and  loaded 
by  transverse  uniform  circular  loadings 
of  arbitrary  radius  and  magnitude.  A 
superposition  routine  allows  positioning 
of  loads  to  simulate  aircraft  gear  con- 
figurations, half  tracks,  or  other  heavy 
equipment  vehicles.  The  full  details  of 
the  development  and  limitations  of  the 
code,  as  well  as  a user  manual  and 
example  input,  are  given  in  Appendix  A. 

Of  primary  importance  is  that  the 
finite-element  code  does  not  embody  any 
assumptions  deviating  from  classical 
elasticity  - aside  from  the  fundamental, 
assumed,  displacement  field  over  each 
element,  which  is  characteristic  of  the 
finite-element  technique.  However,  it 
has  been  shown  that,  as  the  number  of 
elements  is  increased,  the  solution  will 
converge  to  the  exact  solution;  there- 
fore, by  successive  trials  an  "exact-as- 
you-wish"  solution  may  be  obtained. 

As  a check  on  the  computer  code 
the  strain  results  of  the  sea-ice  plate 
studies  are  compared  with  analytical 
solutions  from  the  finite-element 
method.  A description  of  the  laboratory 
plate  experiments  has  already  been 
presented  in  a previous  section. 


Consequently,  only  details  pertinent  to 
strain  measurement  will  be  discussed. 
The  four  gages  (Figure  25)  were  grouped 
in  two  pairs,  so  that  each  pair  measured 
strain  in  two  orthogonal  directions  and 
their  average  cancelled  any  anomalies 
arising  from  asymmetric  responses. 
Strains  were  continuously  recorded,  and 
peak  values  were  easily  distinguished. 
This  pattern  was  repeated  for  different 
peak  load  values.  The  top  gages  in  the 
compression  zone  exhibited  linear 
behavior,  whereas  the  tension  gages  sug- 
gested a nonlinear  response,  perhaps  due 
to  tensile  cracking  of  the  sea  ice. 

Four  constant-temperature  tests 
were  conducted  on  a single  4.5-in.  plate 
at  the  temperatures:  -4,-10,-20,  and  -30 
C.  Figure  41  shows  the  average  experi- 
mental radial  strains  for  both  gages 
located  in  the  tension  and  compression 
zones  at  each  test  temperature.  The 
analytical  strain  distributions  predicted 
by  the  finite-element  technique  for  each 
plate  test  are  superimposed  on  each 
figure.  For  the  analytical  approach,  the 
sea  ice  was  assumed  to  be  homogeneous, 
isotropic,  and  elastic,  with  Poisson's 
ratio  taken  as  0.3  for  all  temperatures. 
Young's  modulus  was  determined  from 
plate-theory  equations  by  using  the 
experimental  strains  as  input.  The  exper- 
imental strains  deviated  somewhat  from 
plate-theory  kinematics  in  that  the 
experimental  strain  values  at  the  tension 
and  compression  zone  did  not  form  a 
straight  line  through  the  plate  centroid. 
Therefore,  in  determining  Young's 
modulus,  more  weight  was  given  to  the 
compression  gage  as  it  was  postulated 
that  the  tension  strains  were  influenced 
by  cracking.  The  computed  modulus  is 
noted  in  each  figure. 

Examination  of  the  results  of  the 
constant-temperature  plate  tests  yields 
the  following  observations. 
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Figure  41a.  Strain  profiles  on  centerline  of  4-1/2-jn.  sea-ice  plate  at  a temperature 
of  -4C. 
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Figure  41b.  Strain  profiles  on  centerline  of  4-1/2-in.  sea-ice  plate  at  a temperature 
of  -IOC. 
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Figure  41c.  Strain  profiles  on  centerline  of  4-1/2-in.  sea-ice  plate  at  a temperature 
of-20C. 
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Figure  41d.  Strain  profiles  on  centerline  of  4-1  /2-in.  sea-ice  plate  at  a temperature 
of  -300. 
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1. The  analytical  approach  is  in 
excellent  agreement  throughout  the  ice 
thickness,  even  in  the  top  fiber 
compression  zone.  This  agreement  is  to 
be  expected,  since  the  plate  dimension 
ratio  - diameter-to-thickness  - is  on  the 
order  of  10,  which  is  the  rule-of-thumb  if 
thin-plate  theory  is  to  be  valid  for  plates 
of  finite  extent. 

2.  It  can  be  observed  that  the 
experimental  strains  in  the  tension  zone 
are  consistently  larger  than  the  analyti- 
cal predictions,  v/hich  motivated  the 
contention  that  the  tension  gages  experi- 
enced a nonlinear  response,  perhaps  due 
to  tension  cracks. 

3.  It  is  confirmed  that  the  calcu- 
lated moduli  demonstrate  a consistent 
relationship  with  the  plate  temperature, 
such  that  colder  temperatures  produce 
higher  moduli. 

Two  experiments  were  conducted  on 
sea-ice  plates  with  linear  temperature 
gradients  induced  across  the  plate  thick- 
ness; -1 1 to  -6C  and  -19  to  -9C  from  top 
to  bottom,  respectively.  Otherwise,  the 
tests  were  identical  to  those  with 
constant-temperature  plates.  Figures  42a 
and  b show  the  average  measured  radial 
strains  along  with  the  radial  strain 
profile  from  the  finite-element  solution. 

Again,  in  determining  the  analytical 
solution,  the  sea  ice  was  assumed  iso- 
tropic and  elastic  with  Poisson's  ratio  set 
at  0.3;  however,  the  modulus  was 
assumed  to  vary  linearly  across  the  plate 
thickness.  In  order  to  determine  this 
linear  distribution,  the  unique  modulus 
distribution  which  produced  a strain  pro- 
file capturing  both  the  compressive  and 
tensile  experimental  strains  was  chosen. 
These  modulus  distributions  are  given  in 
their  respective  figures  in  terms  of  the 
maximum  and  minimum  moduli  at  the 
top  and  bottom  of  the  ice  plate. 


Viscoelastic  Finite-Element 
Technique 

Viscoelastic  materials  are  often 
called  "memory"  materials;  that  is,  the 
stress  in  the  material  is  determined  not 
only  by  the  current  state  of  deformation 
but  also  by  all  past  deformation  states. 
Moreover,  the  "memory"  exhibits  a 
fading  phenomenon  in  that  past  deforma- 
tion states  influence  the  current  stress 
state  to  a lesser  degree  than  do  more 
recent  deformation  states. 

It  is  known  that  sea  ice  behaves 
elastically  only  for  small  to  moderate 
load  magnitudes  applied  over  short  time 
periods.  To  expand  the  range  of  load 
magnitudes  and  time  durations,  the 
theory  of  linear  viscoelasticity  (26,  27) 
has  been  coupled  with  the  finite  element 
technique  to  provide  more  realistic 
characterization  of  sea-ice  behavior. 

A general  summary  of  the  viscoelas- 
tic finite  element  formulation  (5)  is 
presented  in  this  section  while  details  of 
the  development,  as  well  as  computer 
code  user  instructions,  appear  in 
Appendix  B.  The  general,  isotropic,  two- 
dimensional  viscoelastic  stress-strain 
relationship  is  presented  as  a matrix  of 
Stieltjes  integrals  along  with  a special 
characterization  of  the  material 
relaxation  functions  known  as  a prony 
series.  Then  the  viscoelastic  constitutive 
relation  is  incorporated  into  a finite- 
element  formulation  for  axisym- 
metric/plane  strain  geometry  (5).  The 
scope  of  the  mechanics  is  limited  to 
isothermal  linear  viscoelasticity,  small- 
deformation  theory,  and  quasi-static 
loading.  However,  the  inclusion  of 
inertial  terms  for  dynamic  analysis  is  a 
straightforward  extension  of  this  work. 
The  principle  of  virtual  work  is  utilized 
to  develop  a displacement  formulation 
for  the  viscoelastic  model  for  a generic 
class  of  two-dimensional  finite  elements. 
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Depth  (in.) 


Strain  (pin./in.  per  unit  load) 

-.06  -.05  -.04  -.03  -.02  -.0»  0.0  .01  0.2  0.3  0.4  0.5  0.6 


Figure  42a.  Strain  profiles  on  centerline  of  4-1/2-in.  sea-ice  plate  at  a temperature 
gradient  of  -11  to  -6C. 


Strain  (pin./in.  per  unit  load) 

■.06  -.05  -.04  -.03  -.02  -.01  0.0  . 01  .02  . 03  . 04  . 05  . 06 


Figure  42b.  Strain  profiles  on  centerline  of  4-1 /2-in.  sea-ice  plate  at  a temperature 
gradient  of  -19  to  -9C. 


45 


A fluid  foundation  is  incorporated  into 
the  formulation  for  plate-type  structural 
systems,  as  well  as  a superposition 

technique  for  handling  arbitrary  surface 
loads  on  layered  half-space  systems. 

The  finite-element  formulation  re- 
quires that  time-dependent  sea-ice 
properties  be  given  in  the  form  of  relax- 
ation functions  for  both  the  bulk  and 
shear  moduli,  K(t)  and  G(t). 

n -t/X; 

K(t)  = K + E K.  e (21) 

° i = l ‘ 

n -t/Bj 

and  G(t)  = + E G.  e (22) 

o 1 

1=1 

Here  the  relaxation  functions  are  repre- 
sented by  a constant  and  a series  of 
exponential  terms  where  K^,  Kj , . . . , 
Kn , Go , Gj  , . . . , G are  positive 
constant  moduli  value  with  units  of 
stress,  and  X|>  X2>***»Xn’Bl>  62’*** 
Bp  are  positive  constants  called  "relax- 
ation times." 

In  this  prony  series  form,  the  relaxa- 
tion functions  are  capable  of  represent- 
ing a large  class  of  linear  viscoelastic 
materials.  In  essence  any  linear 
viscoelastic  material  which  exhibits  an 
instantaneous  elastic  deformation 
followed  by  a creep  phase  can  be 
modeled  by  Equations  21  and  22  with  a 
sufficient  number  of  exponential  terms. 
In  practice,  the  use  of  more  than  four 
exponential  terms  is  seldom  justified; 
usually  two  terms  are  sufficient. 

Relaxation  Functions  From  Creep 
Data.  Strain  histories  from  previously 
discussed  uniaxial  and  beam  experi- 
mental effort  must  be  converted  into 
relaxation  functions.  First,  all  creep  data 
should  be  discretized  into  data  points  at 
specified  times.  Then,  the  objective  is  to 


fit  a creep  function,  3{t),  in  the  form  of 
a prony  series. 


m -t/  X' 

3{t)  = :Iq  + E :ije  ‘ (23) 

where  3(t)  is  defined  as  the  strain  history 
response  to  a prescribed  unit  Heaviside 
stress  input.  This  curve-fitting  technique 
can  be  accomplished  by  a linear  least- 
square  error  analysis  based  on  m + 1 
unknown  creep  moduli,  3o>  3 | , . . . , 3^^, 
with  preselected  retardation  times,  X-  , 
which  can  be  optimized  by  trial  and 
error.  Once  3(t)  is  determined,  using 
LaPlace  transformation  - the  inverse  of 
3(t),  the  relaxation  function,  Y(t),  can  be 
found.  However,  while  Y(t)  is  of  the  form 
of  Equation  23,  it  does  not  represent 
either  the  bulk,  K(t),  or  shear,  G(t), 
relaxation  functions,  but  rather  Y(t)  = 
E(t),  the  Young's  modulus  relaxation 
function.  Assuming  Poisson's  ratio,  p , to 
be  constant,  then  the  relaxation  func- 
tions for  the  bulk  and  shear  moduli  are 
easily  determined  by  the  following 
relations: 


K(t) 


E(t) 

3(1  - 2p) 


G(t) 


E(t) 

2(1  + p) 


(24) 


Ice-sheet  bending  creates  a critical 
biaxial  stress  state  beneath  the  load 
where  the  ice  is  in  tension  applied 
perpendicular  to  normal  crystal  growth 
(horizontally).  Consequently,  relaxation 
functions  were  determined  for  the 
uniaxial  tensile  laboratory  creep  tests  on 
horizontal  specimens  (see  Figure  32). 
Each  ice  temperature  has  its  own 
governing  relaxation  function  for  the 
modulus  of  elasticity,  E(t). 
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-4°C:  E(t)  = 


ice  surface  temperature  of 
becomes 


* (25)  ^ 5.23(,o’)e-'/‘-“ 


■16C 


-10°C:  E(t)  = 


(31) 


+ 1.07(10^)e'^^^^^°  (26) 

-20°C:  E(t)  = 3.79(10^)e'^^^*^* 

+ l.47(10*’)e'^^‘^^°  (27) 

-27°C:  E(t)  = 5.66(10^)e'^^*'‘^ 

+2.26(10^)e'^^^°*°  (28) 

where  the  term  coefficients  are  given  in 
Ib/sq  in.  while  the  exponential  coeffi- 
cients are  "relaxation  times"  in  minutes. 

The  same  format  was  followed  to 
calculate  the  relaxation  function  E(t)  for 
both  -<»  and  -IOC  laboratory  beam  creep 
data  (shown  in  Figure  35) 

-<»°C:  E(t)  = l.71(10^)e'^^^’*° 

H 5.39(10^)e'^^^^‘'°  (29) 

-10°C;  E(t)  = 2.55(10^)e‘^^^*^^ 

> 9.79(l0*')e'^^^^°  (30) 

Again,  the  units  for  both  term  and 
exponential  coefficients  remain  the 
same. 

A composite  creep  function  was 
generated  for  field  beams,  representing 
the  deflection-time  history  of  beams  3 
and  7 from  Figure  39.  The  resulting 
relaxation  function  for  E(t)  at  an  average 


Good  agreement  exists  between  both  sets 
of  laboratory  results  and  the  field  effort. 
In  general,  the  field  beams  responded 
with  higher  relaxation  values  for  equiva- 
lent ice  temperatures,  but  these  were 
averages  not  taking  the  entire  thermal 
gradient  of  the  beam  into  account. 

Viscoelastic  Ice  Plate  Ex,.^mple.  To 

demonstrate  the  capability  of  the  visco- 
elastic finite-element  computer  code  (5), 
the  structural  response  of  an  infinite 
sea-ice  plate  on  a fluid  foundation  is 
investigated  for  a simulated  aircraft 
loading  condition  and,  using  relaxation 
functions,  is  determined  from  experi- 
mental creep  data.  The  time-dependent 
surface  deflections  of  an  Arctic  sea-ice 
sheet  were  measured  during  a CEL  field 
effort  to  simulate  long-term  (<»  to  6 days) 
loading  of  a C-130  aircraft  on  floating 
ice  sheets. 

The  main  landing  gear  were  repre- 
sented by  four  wood  blocks,  each  having 
the  same  area  as  each  tire  print  of  a C- 
130.  These  blocks  were  placed  at  the 
corners  of  a 5 x Kn  ft  rectangle,  the 
same  spacing  as  the  aircraft  wheels.  A 
series  of  900-gallon  (3.4  kiloliter)  fuel 
bladders  were  stacked  on  a deck  spanning 
across  the  wood  blocks  and  filled  with 
seawater  by  a 600-gpm  pump.  The  test 
setup  is  shown  in  Figure  43.  Of  course, 
the  entire  weight  of  a fully  loaded  C-130 
could  not  be  placed  on  the  loading 
platform,  but  a P/h^  ratio  of  3<60  psi  was 
attained  for  the  relatively  thin  ice  sheet, 
ranging  from  20  to  24  in.  This  ratio 
corresponds  to  a C-130,  weighing  180,000 


Figure  43.  Test  setup  for  simulatirtg  a parked  C-130  aircraft  on  an  annual  ice  sheet. 
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lb,  parked  on  a 55-ir. -thick  ice  sheet. 

No  radial  cr  circumferential 
cracking  was  observed,  but  ice  deflec- 
tions were  measured  beneath  the  load  by 
means  of  a manometer  containing 
antifreeze  and  radially  outward  from  the 
load  by  taking  level  readings.  To  deter- 
mine the  actual  applied  load,  all 
platform  components  were  weighed,  and 
a flowmeter  monitored  the  number  of 
gallons  of  seawater  pumped  into  the 
bladders.  Temperature  and  salinity  pro- 
files were  measured  several  times  during 
the  test  period.  Salinities  remained  fairly 
constant,  averaging  between  5.5  to  6.0 
ppt,  while  average  ice-sheet  surface 
temperatures  ranged  from  -8  to  -30C. 

Two  of  the  deflection  time  histories 
from  the  ice-sheet  creep  tests  are  shown 
in  Figure  <»<»  for  two  different  loadings 
and  ice  thicknesses.  The  solid  lines  are 
results  of  the  finite  element  analysis. 


The  thin  ice  sheet  of  Test  2 was  approxi- 
mated by  a single  layer  with  an  average 
temperature  of  -IOC.  The  relaxation 
functions  for  bulk  and  shear  moduli  were 
calculated  by  substituting  Equation  26 
into  Equation  24.  In  Test  3 the  ice  sheet 
had  a significant  temperature  gradient 
across  the  thickness,  so  two  layers  were 
used  in  the  analysis;  the  top  layer  having 
a mean  temperature  of  -24C  and  the 
bottom  -IOC.  To  model  the  time-depen- 
dent material  properties  of  this  ice 
sheet,  two  sets  of  bulk  and  shear  relaxa- 
tion functions  were  computed  from 
Equations  26  and  28  obtained  from 
uniaxial  laboratory  studies.  Both 
experimental  and  analytical  results 
represent  v/heel  deflections,  including 
the  effect'  of  superposition  by  the  other 
three  wleels.  Both  analyses  agree  with 
field  measurements  very  closely,  not 
differing  by  more  than  0.5  in.  during  80 
hours  of  testing. 
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bearing  capacity  of 
floating  ice  sheets 


Failure  Criteria 


Once  field  and  laboratory  exocri- 

mentation  establishes  reasonable  vaTues 

*or  sea-ice  mechanical  properties  and 

analytical  approaches  have  been  devised 

to  adequately  model  sea-ice  behavior  it 

■s  possible  ,0  solve  real  ice-sl^' 

problems,  such  as  landing  and  paS 

aircralt,  moving  vehicular  trail,?,  and 

fi^rgo.  A rational  utilization  ol 

« ,ng  ,ce  plates  lor  all  these  opera- 

bMHno''’''''”  Iftwwledge  ol  their 

star  an?,?"”''  '"counter 

Start-  and  long-term  loading.  In  the  oast 

mimerous  attempts  have  taen  .maT  to 

Jtermine  the  bearing  capacity  of  ice 

sheets  subjected  to  vertical  loads.  A 

has*l^  '“erature  survey  of  the  subject 
has  been  written  by  Kerr  (28).  ’ 


^'earing  capacity  of  a 
floating  ice  sheet  subjected  to  short- 

sST  V*  ^ maximum 

determined 

given  lead  and  assuming  that  the 
‘t  Next,  the  load,  P 

at  which  the  first  crack  occurs  is  found 
using  the  criterion 


°max  ■ ° f 


(32) 


where  oj  is  generally  given  as  an  exoeri- 
mental  tensile  or  flexural  failure  st?^ss. 

plete  breakthrough  failure  of  the  ice 

fS’t  in  ^'er  tloes  provide 

built-in  factor  of  safety.  Usually, 


<»9 





however,  an  allowable  load,  P^ji,  is 
calculated  as  some  percentage  of  Pj-r  to 
provide  an  additional  safety  factor. 

To  support  loads  for  extended 
periods  of  time,  the  theory  of  elasticity 
is  no  longer  valid  for  analysis.  Applying 
the  theory  of  linear  viscoelasticity, 
deflection-time  histories  can  be  predic- 
ted. One  common  criterion  states  the 
maximum  deflection  should  not  exceed 
the  freeboard  of  the  ice  sheet,  generally 
taken  as  0.08  h,  where  h is  the  ice-plate 
thickness  (29).  If  the  freeboard  is 
exceeded,  water  can  flood  the  surface 
through  cracks  and  disrupt  activities. 
This  criterion  is  definitely  conservative 
and  safe,  and  a maximum  deflection 
criterion  of 


is  used  to  determine  allowable  parking 
times  for  vehicles  and  aircraft  on 
floating  ice  sheets.  Even  though  Equation 
33  represents  a 25%  increase  over  the 
freeboard  criterion,  Kingery  (30)  has 
indicated  there  was  no  failure  during 
field  tests,  even  for  deflection  ratios, 
w/h,  greater  than  0.25. 

Computer  Code  Input 
Parameters 

Elastic  Finite  Element  Method.  The 

input  requirements  for  the  elastic  finite- 
element  code  are  based  on  a circular, 
plate-shaped  solid,  supported  by  a fluid 
foundation.  The  ice  sheet  can  be  of 
arbitrary  dimensions  with  as  many  as  12 
material  layers  and  can  accommodate 
circular  loadings  of  uniform  pressure  but 
arbitrary  radius  and  magnitude.  These 
input  parameters  can  be  divided  into 
three  basic  categories;  (1)  program  con- 
trol and  finite-element  mesh  generation, 
(2)  ice-sheet  material  properties,  and  (3) 
aircraft  or  vehicle  loading  conditions. 


The  first  category  sets  the  boundary 
conditions  and  defines  the  number  of 
material  layers  and  number  of  element 
rows  for  each  layer  of  the  ice  sheet. 
Also,  the  mesh  generator  establishes  the 
number  of  total  elements,  elements 
under  the  load,  and  elements  in  contact 
with  the  supporting  fluid. 

In  the  second  category,  temperatures 
are  assigned  to  the  top  and  bottom  of 
different  material  layers  within  the  ice 
sheet,  along  with  a corresponding  elastic 
modulus  and  Poisson's  ratio  for  each 
temperature.  The  thickness  for  each 
material  layer  is  given,  as  well  as  the 
density  of  the  fluid  foundation  (in  this 
case,  seawater). 

For  the  last  category,  the  input 
parameters  are  the  number  of  circular 
loadings  (e.g.,  wheels  or  tracks)  and  the 
total  load  (weight  of  aircraft  or 
vehicles),  in  addition  to  such  individual 
wheel  or  track  circle  characteristics  as 
(1)  the  radius  and  pressure  of  load;  (2) 
the  x-y  coordinates  (to  locate  each 
loading);  and  (3)  the  fraction  of  total 
weight  carried  by  each  individual 
loading.  Details  of  all  of  these  input 
parameters  are  presented  in  Appendix  A. 

Viscoelastic  Finite  Element  Method. 

The  input  requirements  for  the  visco- 
elastic computer  code  are  similar  to 
those  for  the  elastic  code,  except  for  the 
consideration  of  time  dependency  (5). 
Again,  the  infinite  ice  sheet  is  assumed 
to  be  a large  circular,  plate-shaped  solid 
of  arbitrary  dimensions  resting  on  a fluid 
foundation.  The  ice  sheet  can  have  as 
many  as  12  layers  with  separate  material 
properties  and  be  subjected  to  any 
number  of  circular  loadings  of  uniform 
pressure,  but  with  arbitrary  radius  and 
magnitude.  The  input  parameters  fall 
into  five  different  categories:  (1) 
program  control  and  finite  element  mesh 
generation,  (2)  time-dependent  material 
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Table  1.  Annual  Ice-Sheet  Mechanical  Properties 


Thermal 

Period 

Surface 

Temperature 

("O 

Elastic  Modulus 
Gradient 
(psi  X 10^) 

Flexural 

Strength 

(psi) 

1 

-20  to -10 

7.0  to  2.9" 

70.0 

3 

-10  to  -5 

6.0  to  2.9" 

62.0 

3 

-5  to  -3 

4.0  to  2.9" 

S8.0 

4 

-3  to  -2 

2.9" 

40.0 

"corresponds  to  temperature  gradient,  top  to  bottom. 
^Isothermal  condition. 


properties  for  each  ice-sheet  layer,  (3) 
time  information,  (4)  aircraft  or  vehicle 
loading  conditions,  and  (5)  output  control 
and  suppression. 

As  before,  the  first  category 
established  the  boundary  conditions  and 
defines  the  number  of  material  layers 
and  element  rows  for  each  layer  of  the 
ice  sheet.  The  mesh  generator  sets  the 
number  of  total  elements,  elements 
under  the  load,  and  elements  in  contact 
with  the  fluid  foundation. 

in  the  second  category,  modulus 
coefficients  and  relaxation  times  are 
input  for  both  exponential  series  defining 
the  bulk  and  shear  moduli  for  each  ice 
layer.  The  thickness  of  each  layer  is 
given  as  well  as  the  density  of  the 
supporting  fluid  foundation. 

In  the  tliird  category,  both  the  time- 
step  increment  and  total  number  of  time 
steps  must  be  input,  as  well  as  the  load 
duration. 

For  the  fourth  group  the  input 
parameters  arc  the  number  of  circular 
loadings,  total  load,  and  individual  tire- 
print  characteristics;  (1)  the  radius  and 
pressure,  (2)  tne  x-y  coordinates  (to 
locate  each  loading  on  the  ice  sheet),  and 
(3)  its  percentage  of  the  total  weight. 

The  fifth  category  simply  gives  the 
user  the  opportunity  to  control  the 
printed  output  (deflections,  stresses,  and 
strains).  All  of  these  input  parameters 


are  explained  in  more  detail  in  the  user 
instructions  section  of  Appendix  B. 

Aircraft  Load  Curves 

Calculating  the  required  effective 
ice-sheet  thicknesses  for  the  load  ranges 
of  different  aircraft  requires  exercising 
the  elastic  finite  element  computer  code 
in  two  series  of  repeated  steps  (31). 
Initially,  the  input  parameters  for  a 
single  thermal  period  are  selected  from 
Table  1,  and  several  ice-sheet  thick- 
nesses are  chosen  for  this  period.  By 
performing  consecutive  computer  runs, 
maximum  tensile  stresses  are  found  at 
the  bottom  of  the  ice  sheet  for  each 
specified  thickness.  This  iterative 
process  generates  an  aircraft-load  versus 
ice-thickness  curve  for  a given  thermal 
period,  and  the  entire  cycle  is  repeated 
for  the  remaining  periods.  In  Figure  45 
curves  for  both  C-130  and  C-141  mili- 
tary aircraft  are  shown  for  all  four 
thermal  periods  defined  in  Table  1. 
Period  1 curves  for  both  aircraft  are 
determined  by  comparison  of  the 
calculated  stresses  with  an  allowable 
stress  that  is  found  by  reducing  beam 
failure  strengths  by  30%  (an  equivalent 
safety  factor  of  1.20).  All  other  curves 
reflect  a reduction  of  25%  in  the  failure 
strength,  representing  a safety  factor  of 
1.15. 

Vehicle  Ice  Thickness  Table 

Minimum  ice-sheet  thicknesses  for 
typical  polar  vehicles  are  presented  in 
Table  2 (31).  Again,  the  iterative  process 
is  repeated  for  each  thermal  period. 
However,  for  relatively  constant  loads 
(like  vehicles),  the  iterations  converge  to 
a single,  tabular  point  rather  than 
forming  a curve  as  in  the  case  of 
aircraft,  which  have  a wide  cargo/fuel 
range.  All  vehicles  are  a sampling  of 
typical  logistics  support  equipment 
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operating  in  polar  regions.  Vehicle 
weights  used  in  thickness  calculations 
are  also  given  in  Table  2.  Allowable 
stresses  represent  a sufficient  reduction 
in  failure  strength  to  produce  a safety 
factor  of  1.50. 

Allowable  Parking  Times  for 
Aircraft 

Maximum  allowable  parking  times  in 
hours  for  both  C-130  and  C-141  aircraft 
are  presented  in  Table  3.  Using  the 
viscoelastic  finite  element  computer 
code  (5),  time-dependent  deflections 
beneath  the  main  landing  gear  of  both 
aircraft  were  calculated  for  three 
separate  ice-sheet  thicknesses  at  each 
thermal  period  of  operation.  The  first  ice 
thickness  listed  for  each  period  was 
taken  from  the  minimum  ice  thickness 
curves  in  Figure  45  for  the  particular 
aircraft  weight.  An  aircraft  can  hardly 
be  allowed  to  park  on  an  ice  sheet 
thinner  than  landing  thickness.  The  re- 
maining two  ice  thicknesses  were  simply 
6 and  12  in.  greater  than  the  first.  The 
four  seasonal  periods  are  defined  by  the 
surface  ice  temperatures  listed  in  Table 
I,  but  the  elastic  modulus  gradients 
presented  in  Table  1 arc  invalid  for  these 
time-dependent  cases.  Thus,  for  thermal 
period  1 the  ice  sheet  was  divided  into 
four  layers,  each  governed  by  its  own  set 
of  relaxation  functions  for  bulk  and  shear 
moduli  calculated  from  Equations  25 
through  28.  For  period  2 the  ice  sheet 
was  divided  into  three  layers,  using 
relaxation  functions  based  on  Equations 
25  through  27.  This  pattern  was  repeated 
for  period  3 using  three  layers,  while  the 
ice  sheet  was  considered  virtually  iso- 
thermal for  period  4,  governed  only  by 
Equation  25.  The  failure  criterion  in 
Eqijation  33  was  applied  to  determine  the 
allowable  parking  times.  Superposition 
was  employed  to  consider  contributions 
from  all  other  aircraft  wheels.  It  should 


be  noted  that  during  these  computer 
simulations  of  parked  aircraft,  the 
principal  biaxial  tensile  stresses  at  the 
bottom  of  the  ice  sheet  relaxed  on  the 
order  of  60%  of  their  initial  value  at  t = 
80  hr.  The  fact  that  stresses  are  reduced 
as  time  goes  on  should  encourage  the 
development  of  a failure  criterion  for 
sea-ice  sheets  based  on  limiting  strain. 

Operational  Field  Procedure 

Once  ice-thickness  curves  and  tables 
are  determined,  proper  field  procedure  is 
still  necessary  to  operate  safely  on  any 
ice  sheet.  Of  course,  both  ice  thicknesses 
and  temperatures  must  be  monitored.  At 
the  start  of  the  operating  season,  a 
detailed  thickness  survey  should  be  made 
down  each  edge  of  sea-ice  runways, 
alongside  ice  roads,  and  around  all  other 
operating  areas.  Maintenance  of  a com- 
plete up-to-date  temperature  record  of 
the  ice  sheet  is  necessary  for  a verifica- 
tion that  the  proper  thermal  period  is 
being  used.  Ice  temperatures  should  be 
measured  by  thermocouple  or  thermistor 
stations  at  five  or  more  locations 
dispersed  for  general  coverage  of  the 
entire  operational  area.  Routine 
inspection  of  the  ice-sheet  surface 
should  be  made  for  existence  of  cracks. 
Such  cracks  should  be  recorded  for 
length,  location,  and  whether  they  are 
wet  or  refrozen. 


CONCLUSIONS 

1.  The  comprehensive  testing  program  on 
the  mechanical  properties  of  sea  ice  has 
provided  a wealth  of  basic  information 
necessary  for  understanding  sea-ice 
behavior  under  different  loading  condi- 
tions. 
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Table  2.  Minimum  Allowable  Ice  Thicknesses  for 
Representative  Equipment  and  Vehicles 


Vehicle 

Minimum  Allowable  Ice  Thickness  (in.) 

or 

Equipment 

Thermal 

Thermal 

Thermal 

Thermal 

Period  1 

Period  2 

Period  3 

Period  4 

1.  Pickup,  IT-Dodge  W300 
(GVW  = 9,000  lb)** 

11 

13 

17 

21 

2.  Cargo/personnel-Nodwell 

45TT  (15,000  lb) 

16 

18 

23 

28 

3.  Personnel  carrier-Nodwell 

17 

21 

28 

34 

RNllO  (GVW  = 27,500  lb) 

4.  Personnel  carrirr-Trackmaster 

10 

12 

16 

20 

601  (7,540  lb) 

5.  Caterpillar  D-4  STD 
(17,250  lb) 

20 

24 

31 

38 

6.  Caterpillar  D-8  l.GP 
(85.0001b) 

44 

48 

59 

71 

7.  C^aterpillar  950 
(24.900  lb) 

24 

28 

35 

42 

8.  Caterpillar  955  STD 
(30.200  lb) 

79 

33 

42 

50 

9.  Road  grader  - Cat  12  F 
(28.310  lb) 

26 

29 

37 

44 

10.  Crane.  wheeled-Pettibone  70 
(70.7301b) 

40 

45 

57 

69 

11.  5T  tractor  (military)  ♦ I2T 

trailer  (FruehauO 
(GVW  = 56.0001b) 

27 

31 

39 

43 

12.  lOT  tractor  (Fabco)  ♦ 20T 

trailer  (Fabco) 

(GVW  * 65.0001b) 

29 

33 

42 

51 

13.  lOT  tractor  (military)  ♦ 60T 

lowboy  (military) 

(GVW  = 180.000  lb) 

52 

58 

74 

88 

14.  Caterpillar  D-8  l.GP  ♦ 20T 
sled  (GVW  = 145.000  lb) 

45 

51 

61 

72 

‘'(iVW  is  gross  vehicle  wvight. 
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Table  3.  Allowable  Parking  Times  for  C-1 30  and  C-141  Aircraft  for  Different  Ice-Sheet  Thicknesses 


Thermal 

Period" 


Icc-Sheet 

Allowable 

Icc-Shcct 

Allowable 

Thickness 

Parking  Times 

Thickness 

Parking  Times 

(in.) 

(hr) 

(in.) 

(hr) 

C-1  30  Aircraft'’ 

C-141  Aircraft'" 

54 

7 

74 

15 

60 

16 

80 

31 

66 

37 

86 

43 

60 

8 

78 

14 

66 

21 

84 

28 

72 

42 

90 

44 

76 

55 

102 

106 

82 

87 

108 

121 

88 

122 

114 

152 

90 

97 

122 

149 

96 

130 

128 

180 

102 

164 

134 

208 

"Defined  by  ice  temperatures  from  Table  1. 
'’Aircraft  weight.  180.000  Ib. 

“^Aircraft  weight.  300.000  Ib. 


2.  Two  finite-element  computer  codes 
developed  at  CEL  have  established  a 
capability  to  analyze  ice-sheet  behavior 
through  the  v.scoelastic  range.  This 
analytical  software  has  the  ability  to 
consider  the  influence  of  the  existing 
temperature  gradient  across  the  ice 
sheet  and  to  handle  multiple  wheel  or 
track  loadings  on  the  ice  sheet. 

3.  By  a combination  of  mechanical 
property  testing  and  finite-element 
analysis,  workable,  realistic  ice  thickness 
requirements  have  been  determined  for 
aircraft  and  vehicle  operations  on 
floating  sea-ice  sheets,  as  well  as 


allowable  times  for  parked  aircraft  and 
effects  of  cargo/fuel  storage  and  snow 
overburden  near  ice-sheet  operations. 


RECOMMENDATIONS 

1.  Further  experimental  and  analytical 
effort  is  recommended  to  develop  an 
accurate,  sensible  failure  criterion  to 
better  define  the  safe  limits  of  long- 
term sea-ice-sheet  behavior. 

2.  Any  further  laboratory  or  field  testing 
of  sea  ice  will  require  better  strain  and 
deflection  measurement  techniques  to 
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insure  accurate  and  reliable  deformation 
data  necessary  to  establish  any  failure 
criterion  or  to  verify  developed  constitu- 
tive relations. 

3.  It  is  recommended  that  a user  manual 
be  prepared  containing  a complete  and 
detailed  description  of  the  application 
and  limitation  3 of  the  finite-element 
computer  codes,  including  a complete 
listing  of  computer  programs,  input 
format,  and  varied  examples  to  illustrate 
the  solution  of  elastic  and  viscoelastic 
ice  problems.  The  manual  preparation 
should  be  preceded  by  a thorough  review 
of  the  computer  codes  for  restructuring 
the  program  toward  simplification  to 
assist  the  practitioner  who  may  not  be 
highly  skilled  in  computer  technology. 

k.  As  an  ultimate  goal,  it  is  recommen- 
ded that  viscoplastic  theory  be  incorpor- 
ated into  a finite-element  program  for 
analyzing  systems  composed  of  sea  ice. 
Viscoplasticity  is  a logical  extension  of 
linear  viscoelastic  theory  such  that 
plastic  flow  is  also  included  in  the 
formulation.  For  example,  viscoplastic 
models  are  capable  of  replicating 
primary,  secondary,  and  tertiary 
deformation  characteristics  of  sea  ice. 
Moreover,  nonlinear  and  time-dependent 
failure  criteria  are  inherent  in  the 
formulation.  Further,  these  failure 
criteria  (i.e.,  flow  laws)  can  be  rationally 
defined  to  exhibit  different  behavior 
characteristics  for  tension  and  compres- 
sion stress  states  as  is  observed  with  sea 
ice. 

Preliminary  investigations  of  apply- 
ing viscoplastic  models  to  sea  ice  have 
been  extraordinarily  encouraging  (32).  It 
is  recommended  that  these  models  be 
vigorously  researched  and  developed  as 
the  final  answer  to  the  characterization 
of  sea  ice. 
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Appendix  A 


FINITE-ELEMENT  COMPUTER  CODE  FOR  ELASTIC  ANALYSIS  [4] 

ANALYTICAL  DEVELOPMENT 
General 


The  main  portion  of  the  general  axisymmetric  finite  element  computer 
code  employed  in  this  study  was  written  by  Edward  L.  Wilson  [25]  and  extended 
by  M,  G.  Katona  [4],  Formulation  is  based  on  the  direct  stiffness  method  with 
triangular  or  quadrilateral  conical  elements  representing  the  axisymmetric  structure. 
The  structural  system  may  be  of  complex  configuration  and  composite  construc- 
tion, and  the  constitutive  relations  may  be  orthotropic,  nonlinear,  and  temperature 
dependent.  External  loading  may  he  mechanical  or  thermal  in  nature  with  an 
option  to  include  body  forces. 

A special  mesh  generator  was  added  to  the  program  to  facilitate  input 
for  the  plate-shaped  solids  considered  in  this  study.  In  addition,  subroutines 
to  simulate  a fluid  foundation  and  superposition  of  loads  from  arbitrary  load 
configurations  were  incorporated  into  the  main  computer  deck.  These  sub 
routines,  as  well  as  input  instructions  and  examples,  are  discussed  in  the 
following  sections. 

Fluid  Foundation 

Development  of  the  simulated  fluid  foundation  was  achieved  by 
calculating  a vertical  resistance  proportional  to  the  displaced  volume  and 
adding  it  to  the  global  stiffness  matrix  of  the  solid  for  every  node  on  the 
solid-fluid  interface.  In  formulating  the  vertical  resistance,  compatibility  of 
displacements  requires  off-diagonal  terms  as  well  as  diagonal  terms,  all  of  these 
terms  are  determined  by  energy  principles  consistent  with  the  finite  element 
formulation.  The  following  derivation  is  based  on  these  concepts. 

Consider  a solid  supported  by  a fluid  foundation  with  a typical  node 
trio  i,  j,  k and  an  arbitrary  vertical  displacement,  Aj,  at  node  j as  shown  in 
Figure  A-1 . Denoting  the  horizontal  surface  of  a one-radian  segment  bounded 
by  r^  and  fj  as  face  A,  and  the  horizontal  surface  bounded  by  fj  and  as  face 
B the  following  relations  are  valid  after  point  j detorms  through  the  imposed 
displacement  Aj. 

Vertical  displacements  consistent  with  the  finite  element  interpolation 
lunction  along  faces  A and  B are  given  by: 
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''A(r)  = fry  Aj  . Ti  < r < Tj 

i '1 

(A-1) 

•"k  - r 

Vb('’)  = ; Aj  , Tj  < r < r,. 

'k  'i 


Vertical  loading  intensities  along  faces  A and  B are  expressed  as; 


1 

fA(r)  = J = Tw^Airlr  , Pj  < r < Tj 

0 

(A-2) 

1 

= j = Tw^etrlr  , Pj  < r < r,, 

0 


Figure  A-1.  Illustration  of  deformation  from  displacement,  Aj. 
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The  net  vertical  restoring  forces  on  faces  A and  B are; 


'"a  = / (4rf  - 2rjri  - 2r2) 

h 

(A-3) 

Fb  = / Vr)dr  = (2r2  + 2r,.rj  - 4r?> 

Next,  energy  principles  are  employed  to  determine  kinematically 
equivalent  nodal  forces,  F^j,  F^j,  F^j  and  Fg,,,  depicted  in  Figure  A-2  with 
the  stipulation  that: 


+ Faj  = F^ 

(A-4) 

Fbj  + Fb,  = Fb 


Figurt  A-2.  Kintmatically  tquivaitnt  nodal  forces. 


Thus,  for  face  A an  impressed  force,  F^j,  is  sought  which  displaces 
a unit  virtual  load  at  node  j through  the  distance  Aj.  The  magnitude  of  this 
force  can  be  determined  by  equating  the  virtual  work  of  the  unit  load  to  the 
work  of  the  virtual  loading  intensity  along  face  A as  it  deforms  under 
impressed  load  The  virtual  loading  intensity  is  simply  the  ratio  I^F^j 
times  the  loading  intensity  given  by  Equation  A-2.  Therefore,  the  virtual 
work  relation  is  given  as: 


V„(r) 


i *a'** 

= f f -r 

I -I 


which  gives; 


Fai  * (3rf  - 2rjri  - rf) 


(A-5) 


By  the  same  procedure  the  kinematically  equivalent  nodal  force 
can  be  determined.  By  applying  Equations  A-3  and  A-4  the  complete  set 
of  nodal  forces  can  be  determined  as: 


C « (.2  -2) 


^ 13-f  - 2-,-,  - -i 


^Bi  “ + 2r,rk  + rj) 


Bi 


12 


(r^  - rf) 


(A-6) 


Equations  A-6  form  the  basis  of  the  algorithm  in  the  computer  code 
for  determining  the  stiffness  additions  representing  the  fluid  foundation. 

Superposition 

The  purpose  of  the  superposition  development  is  to  determine  the 
combined  state  of  stress  in  a plate-solid  resulting  from  an  arbitrary  loading 
contiguraiion  of  circular  loads.  As  is  generally  the  case,  the  prirtciple  of 
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superposition  has  inherent  limitations;  first,  the  constitutive  relations  must 
be  linear,  and  second,  the  plate-solid  must  be  free  of  variations  of  material 
property  and  layer  geometry  with  respect  to  any  horizontal  plane.  Fortu- 
nately typical  ice  sheets  are  composed  of  horizontal  layers  which  have  material 
property  variations  and  terriperature  changes  solely  in  the  vertical  direction. 
Moreover,  the  "infinite"  extent  of  the  ice  sheet  I'ullifies  the  effect  of  bound- 
aries confounding  the  superposition  principle;  hence,  the  only  significant 
limitation  is  thi.  restriction  to  linear  stress-strain  laws. 

The  process  of  determining  the  superimposed  stress  state  begins  by 
choosing  a point  in  the  plate-solid  where  the  combined  state  of  stress  is  to  be 
calculated.  The  chosen  point  is  referenced  by  an  XYZ  rectangular  coordinate 
system  at  (0,  0,  z*)  where  z*  is  an  explicit  depth  from  the  coordinate  origin 
on  the  plate  surface  to  the  designated  point  in  the  plate-solid.  Next,  for  each 
load,  i,  the  load  center  is  referenced  by  the  coordinate  pair  (Xj,  yj),  then 
transformed  to  a local  polar  coordinate  pair  (fj,  Oj)  by  the  transformations: 


d,  = arctan(y,/Xj) 


The  ordered  pair  (rj,  z*)  defines  a corresponding  coordinate  point  in  the 
plate-solid  of  a single-load  finite  element  program  where  stresses  can  be 
extracted  from  data  storage.  However,  the  local  state  of  stress  tfius  obtained 
is  at  angle  0^  with  respect  to  the  global  XYZ  coordinate  system;  consequently, 
a coordinate  transformation  is  performed  to  determine  the  local  state  of  stress 
in  the  global  system.  This  process  is  repeated  and  the  stress  contribution  of 
each  load  at  the  designated  point  is  summed  to  yield  the  complete  combined 
state  of  stress.  Lastly,  the  principal  stresses  and  unit  normals  are  calculated 
from  the  complete  combined  state  of  stress. 


DISCUSSION  OF  APPLICATION 
General 


Due  to  the  general  nature  of  the  finite  element  code  developed  in 
this  study,  a variety  of  engineering  investigations  could  be  undertaken  such 
as  determining  the  effect  of  material  and  geometrical  variations  on  the  struc- 
tural behavior  of  ice  sheets.  However,  the  objective  of  this  Appendix  is  aimed 
at  development  of  t!ie  analytical  model  rather  than  application.  Thus,  the 
limited  number  of  test  problems  solved  in  this  study  are  directed  toward 
representation  of  an  inifinite  boundary  by  a finite  model. 
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Representation  of  an  Infinite  Boundary 


For  many  applications  an  infinite  expanse  of  ice  is  characteristic  of 
the  problem  to  be  solved.  In  such  cases  defining  an  infinite  expanse  by  an 
axisymmetric  solid  of  finite  extent  is  justifiable  only  if  the  boundary  radius 
is  chosen  large  enough  so  that  changes  in  the  boundary  condition  do  not 
appreciably  alter  the  structural  response. 

Preliminary  investigations  show  that  the  structural  response  of  typical 
ice  sheets  tends  toward  plate-theory  behavior  in  the  outer  region  beginning 
several  loading  radii  away  from  the  load  center.  Consequently,  the  minimum 
boundary  ladius  requirv?d  may  be  determined  in  terms  of  the  plate- theory 
parameter,  radius  of  relative  stiffness,  defined  as: 


C 


Eh^ 

12(1 


where  C 

E 

h 

k 


radius  of  relative  stiffness 
Young's  modulus 
thickness  of  plate 
Poisson's  ratio 
fluid  density 


Figure  A-3  illustrates  the  variation  in  maximum  tensile  stress  for  a 
typical  ice  sheet  and  loading  as  a function  of  boundary  length  in  units  of  £. 
The  three  boundary  conditions  specified  on  the  periphery  which  demonstrate 
convergence  are  (1)  completely  free,  (2)  horizontal  fixed  with  vertical  free, 

(3)  completely  fixed.  As  illustrated  by  the  graphs,  convergence  occurs  at 
approximately  1C2  with  the  fixed-boundary  condition  converging  fastest. 
Moreover,  a maximum  lower  and  a minimum  upper  bound  is  given  by  the 
fixed  and  free  boundary  conditions,  respectively.  Thus,  in  general,  a bound- 
ary radius  of  106  with  a fixed  boundary  condition  will  adequately  represent 
an  infinite  ice  sheet  and,  if  desired,  the  accuracy  may  be  tested  by  applying 
a free  boundary  condition  to  bracket  the  solution. 

Although  the  method  presented  above  will  suffice  for  most  applica- 
tions, optimum  use  of  the  element  is  not  being  achieved  in  the  outer  regions 
where  plate  behavior  is  dominant.  A higher  degree  of  accuracy  and  saving  of 
computer  time  could  be  realized  by  combining  the  Ritz  method  with  the 
finite  element  technique  in  order  to  simulate  the  infinite  boundary. 
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Boundary  Radius,  r^,,  in  Units  of  ( 

Figure  A-3.  Variation  of  maximum  tensile  stress  for  increasing  boundary  radius  for 
typical  ice  sheet. 

PROGRAM  INPUT  FOR  PLATES  ON  FLUIDS 

The  input  requirements  are  given  for  a circular,  plate-shaped  solid  sup- 
ported by  a fluid  foundation.  The  plate  solid  may  be  of  arbitrary  dimensions 
with  as  many  as  twelve  material  layers  and  can  accommodate  circular  loadings 
of  uniform  pressure  of  arbitrary  radius  and  magnitude.  A superposition  routine 
is  available  to  determine  the  local  state  of  stress  resulting  from  any  configuration 
of  circular  loadings  at  any  desired  point. 

A.  Identification  Card  {72H) 

Columns  1 - 6 MESHGN  (This  activates  mesh  generator) 

7 — 72  Descriptive  title  of  problem 
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B.  Fluid  Property  Card  (F10.0) 

Columns  1-10  Density  ot  supporting  fluid 


C.  Control  Card  (lOX,  15,  25X,  FI 0.0, 15,  5X,  15) 


Columns  1 1 
41 
51 
61 


15  Number  of  different  material  layers 
(maximum  of  12) 

50  Reference  temperature  (s+ress-free 
temperature) 

55  Number  of  approximations  (set  to  1 for 
elastic) 

65  Boundary  code;  0,  1, 2,  or  3 

0 — Nodes  on  periphery  are  free 

1 — Nodes  on  periphery  are  vertically  free 

2 — Nodes  on  periphery  are  horizontally 

free 

3 — Nodes  on  periphery  are  fixed 


D.  Material  Layer  Properties  — (one  card  per  layer,  followed  by 

temperature  cards) 


-First  Card  (215,  2F10.0, 110,  3F10.0) 


Columns  1 - 5 
6-  10 


11 

- 20 

21 

-30 

31 

-40 

41 

- 50 

51 

-60 

61 

- 70 

Layer  identification  number  (from  1 to  12) 
Number  of  different  temperatures  for  which 
layer  properties  are  given  (8  maximum, 

1 minimum) 

Mass  density  of  material  (leave  blank  when 
not  considering  body  forces) 

Ratio  of  plastic  modulus  to  elastic  modulus 
(leave  blank  for  elastic  problem) 

Number  of  element  rows  in  this  layer 
Thickness  of  this  layer 
Temperature  at  top  of  layer 
Temperature  at  bottom  of  layer 
(Note:  program  will  linearly  distribute 
temperature  through  the  layer) 


-Following  Cards  (one  card-group  for  each  temperature  indicated 
above)  (8F10.0) 


21-30  Poisson's  ratio,  v„ 

31  — 40  Modulus  of  elasticity,  Eg 
41-50  PoiSson's  ratio,  Vg,  and  Vg^ 

51  — 60  Coefficient  of  thermal  expansion,  and  o, 
61  - 70  Coefficient  of  thermal  expansion,  Og 
71-80  Yield  stress  (leave  blank  for  elastic  case) 


E.  Mesh  Generation  Control  Cards  (2F  10.0,  415) 


F. 


Columns  1 — 10 
11-20 
21  -25 
26  - 30 
31-35 


36-40 


Radius  of  loading 
Pressure  of  loading 

Number  of  radial  elements  under  the  load 
Number  of  radial  elements  total 
Number  of  loading  radii  to  periphery  (to 
represent  an  infinite  plate  this  number 
should  be  at  least  ten  times  the  radius  of 
relative  stiffness  divided  by  the  radius  of 
loading) 

Number  of  elements  in  contact  with  fluid 
foi'  'dation  (in  general  set  this  equal  to  total 
number  of  radial  elements) 


Superposition  Control  Cards  (15,  2F10.0) 


Columns  1 - 5 

6-15 

16-25 


Number  of  circular  loadings  (number  of 
wheels) 

Total  load  (weight  of  aircraft) 

Z-ordinateat  which  superposition  is  desired 
(Note:  The  Z-ordinate  in  conjunction  with 
the  arbitrary  surface  coordinates, 

(x,  y)  define  the  spatial  point, 

(x,  y,  z)**(0,  0,  Z-ord.),  at  which 
superposition  will  occur.  The  com- 
plete local  state  of  stress  will  be 
determined  and  principal  stresses 
with  unit  normals  will  be  computed.) 


-Individual  wheel  load  characteristics — one  card  for  each  wheel 


(3F10.0) 

Columns  1 — 10  X-ordinate  of  load  center 
11-20  Y-ordinate  of  load  center 
21  — 30  Fraction  of  total  load  taken  by  this  wheel 
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SAMPLE  INPUT  PROBLEMS 


I .  A single  wheel  load  of  8.1 5 inch  radius  and  200  psi  pressure  acting  on  a 
60-inch-thick  ice  sheet  of  "infinite"  extent. 

Column#!  5 10  15  20  25  30  35  40  45  50  55  60  65 

Card  # 

1 MESHGN  WHEEL  OF  A C-141  PLANE  ON  60  INCHES  OF  ICE 


2 

0.037 

3 

1 

4 1 

1 

10  60.0 

5 

500000. 

0.3 

500000.  0.3 

6 

8.15  200.0 

5 

20  880  20 

7 

1 21600.  -60. 

8 

0.0  0.0 

1.0 

Comments: 

1 . In  card  #3  column  #65,  the  boundary  code  of  "3"  signifies  a fixed 
boundary  which  best  represents  the  infinite  boundary  condition. 

2.  In  card  #6,  the  infinite  boundary  is  represented  by  880  loading  radii 
which  was  determined  from  10  times  the  radius  of  relative  stiffness  divided 
by  the  loading  radius.  ( 1 0X720/8. 1 5) 

3.  The  superposition  cards  (#7  & 8)  are  trivial  in  this  case  since  only  one 
wheel  is  being  considered.  If  the  problem  were  to  determine  the  maximum 
stress  from  the  C-141  gear-group  of  four  wheels,  the  superposition  of  maximum 
tensile  stresses  would  be  specified  as  follows: 

Card  # 


7 

4 

86400. 

-60. 

8 

0.0 

0.0 

0.25 

9 

32.5 

0.0 

0.25 

10 

0.0 

48.0 

0.25 

11 

32.5 

48.0 

0.25 

6S 


II.  A single  C-141  wheel  load  acting  on  a 60-inch-thick  ice  sheet  composed 
of  two  different  layers  of  material  properties  each  30  inches  thick.  There  is 
a temperature  differential  across  each  layer  which  in  turn  affects  the  material 
properties. 


Column  #1 
Card  # 

5 

10  15  20 

25 

30  35  40 

45  50 

55  60 

65 

1 MESHGN 

TWO  LAYERS  OF 

ICE  WITH  TEMPERATURE 

DIFFERENTIAL 

2 

0.037 

3 

2 

1 

3 

4 

1 

3 

5 

30.0 

-20.0 

-10.0 

5 

-25.0 

600000. 

0.3 

600000. 

0.3 

6 

-15.0 

550000. 

0.3 

550000. 

0.3 

7 

-5.0 

520000. 

0.3 

520000. 

0.3 

8 

2 

3 

5 

30.0 

-10. 

-0.0 

9 

-15.0 

510000. 

0.3 

510000. 

0.3 

10 

-5.0 

490000. 

0.3 

490000. 

0.3 

11 

10.0 

460000. 

0.3 

460000. 

0.3 

12 

8.15 

200.0 

5 

20  880  20 

13 

1 21600.  -60. 

14 

0.0 

0.0 

1.0 
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Appendix  B 


VISCOELASTIC  ANALYTICAL  FORMULATION  AND 
COMPUTER  CODE  USER  INSTRUCTION  (5) 


ANALYTICAL  FORMULATION 

The  following  discussion  presents  a brief  description  of  the  finite-element 
formulation  implementing  linear  viscoelasticity  and  its  associated  assumption  and 
limitations,  and  computer  code  user  instructions  developed  by  Katona  (5). 
Viscoelastic  materials  are  often  called  "memory"  materials  because  the  current 
state  of  stress  in  the  material  is  determined  not  only  by  the  current  deformation, 
but  also  by  all  past  deformation  states.  Moreover,  the  "memory"  exhibits  a fading 
phenomenon,  since  past  deformations  have  less  influence  on  the  current  stress  state 
than  do  more  recent  deformation  states. 

Viscoelastic  Representation 

The  above  characteristics  can  be  modeled  either  by  differential  or  integral 
equations  (26).  The  differential  equation  approach  interprets  the  viscoelastic  model 
in  terms  of  a mechanical  analogy  composed  of  springs  and  dashpots.  However, 
integral  equations  are  more  readily  adaptable  to  laboratory  testing  results  and, 
thus,  easily  applied  to  the  solution  of  boundary  value  problems.  An  expression  for 
stress  as  an  integral  function  of  any  arbitrary  strain  input,  e (t),  can  be  written  as: 

o(t)  = Y(t)e^  + r Y(t-T)  ^4^ 

'V 

where  Cq  is  the  initial  elastic  strain,  while  Y(t)  is  termed  the  relaxation  function 
and  is  an  inherent  characteristic  of  the  material  model. 

It  is  desirable  here  to  elaborate  on  the  merits  of  the  integral  form  of  Equation 
B-1.  First,  standard  laboratory  procedures  for  testing  viscoelastic  materials  involve 
only  direct  measurement  of  the  relaxation  function  (or  its  inverse,  the  creep 
function).  A test  specimen  is  subjected  to  a constant  state  of  unit  deformation,  and 
the  resulting  stress  is  measured  as  a function  of  time.  Clearly,  it  is  not  necessary 
to  deduce  an  equivalent  mechanical  analog,  since  the  laboratory  data  can  be  used 
directly  in  Equation  B-1. 

For  isotropic  materials,  it  can  be  shown  that  there  are  only  two  independent 
relaxation  functions  required  to  describe  the  constitutive  relation.  In  solving 
boundary  value  problems,  it  is  convenient  to  choose  the  relaxation  functions  as  the 
responses  to  bulk  and  shear  deformation.  The  bulk  relaxation  function,  K(t),  is 
defined  as  the  hydrostatic  stress  history  due  to  a prescribed  unit  Heaviside  volume 
deformation.  Likewise,  the  shear  relaxation  function,  G(t),  is  defined  as  the  shear 
stress  history  due  to  a prescribed  unit  Heaviside  shear  strain.  Both  K(t)  and  G(t)  are 
positive,  monotonically  decreasing,  independent  relaxation  functions.  They  can  be 
written  in  the  form  of  a prony  series: 
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(B-2) 


ptiq 

K(t)  = Kq  + E G(t)  = G + E G^e'^^ 

i=l  * i-1 

Here,  both  relaxation  functions  are  represented  as  a constant  and  summation  of 
exponential  terms;  where  Kq,  Kj,  . . . , KmK»  Go,  Gi,  . . . , Gm^are  positive 
moduli  constants  with  units  of  stress,  and  X j , • • . > 6j  , . . . >6mG  ^^e 

positive  constants  referred  to  as  "relaxation  times." 

In  the  prony  series  form,  the  relaxation  functions  are  capable  of  representing  a 
large  class  of  linear  viscoelastic  materials.  In  essence,  any  linear  viscoelastic 
material  which  exhibits  an  instantaneous  elastic  deformation  followed  by  a creep 
phase  can  be  modeled  by  Equation  B-2  with  a sufficient  number  of  exponential 
terms.  In  practice,  it  is  seldom  justified  to  use  more  than  four  exponential  terms; 
usually  two  terms  are  sufficient.  It  can  be  noted  that  neither  K(t)  nor  G(t)  is  the 
most  convenient  response  for  laboratory  measurement,  but  converting  normal 
laboratory  results  will  be  discussed  later  for  specific  viscoelastic  saline-ice  tests. 

Displacement  Method 

in  this  derivation  a brief  outline  of  the  finite  element  technique  is  presented, 
wherein  quasi-static  loading,  small  deformations,  axisymmetric/plane  strain 
geometry,  and  the  absence  of  thermal  effects  are  assumed.  First,  the  viscoelastic 
constitutive  relation  of  Equation  B-l  can  be  presented  as  a matrix  of  Stieltjes 
integrals: 


(o(t)|  = ( ^ r|D(t-T)j  (e(t)}dt  (B-3) 

where  {o(t)}  is  the  stress  vector  at  the  current  time,  t;  (e(T)}  is  the  strain  history 
vector  for  0<  t <t;  and  [D(t-T)]  is  the  constitutive  matrix,  composed  of  the  two 
independent  relaxation  functions,  K(t)  andG(t),  in  the  form  of  Equation  B-2. 

When  this  constitutive  model  is  used  in  a displacement  finite-element 
formulation  based  on  virtual  work,  the  resulting  equilibrium  equations  have  the 
form; 

f |K(t-T)|  ^ {u(T)}dTc  (R(t)!  (B-4) 


where  (u(T)}is  the  nodal  point  displacement  history  vector,  {R(t))  is  the  current 
load  vector,  and  [Kft-t)]  is  the  viscoelastic  stiffness-history  matrix.  For  the  sake 
of  clarity,  inertial  terms  are  not  included  in  Equation  Q-k  because  they  can  be 
treated  independently  of  the  viscoelastic  terms. 

In  order  to  solve  Equation  B-«*  for  the  nodal  displacments  at  the  current  time,  t 
::  tp,  a step-by-step  numerical  integration  procedure  must  be  employed.  To  this 
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en  ',  the  time  integral  is  subdivided  into  n time  increments  and  the  time  increment, 
A t,  is  chosen  sufficiently  small  to  assume  that  the  nodal  displacements  vary 
linearly  within  each  time  step.  With  this  approximation,  the  nodal  velocities 
9/  9t{  u(t  )}  can  be  replaced  by  displacement  increments  for,  each  time  step,  At[^, 
by  (1/  At|^)  { Au|^}.  Next,  by  shifting  all  time-interval  integrations  to  the  right-hand 
side,  except  the  current  time  interval.  At,-,,  Equation  can  be  expressed  as 

[K(iit„)]  = {R(t„)}  - (B.5) 

where  [K(Atr,)]  is  the  current  viscoelastic  stiffness  matrix  which  is  dependent 
only  on  the  current  time-step  size,  and  {H(tp|)  } is  the  viscoelastic  force-history 
vector  dependent  upon  all  past  time  steps. 


!»«„)} 


n-1 


I 


(B-6) 


Using  Equation  B-5,  it  is  possible  to  solve  for  { AUp,}  in  a step-by-step  fashion 
beginning  with  elastic  solution,  n = 0,  and  increasing  n successfully  through  the 
time  of  interest  while  accumulating  the  incremental  displacements.  Since  the 
current  viscoelastic  stiffness  matrix  is  dependent  only  on  the  time-step  size,  the 
matrix  need  only  be  reformed  and  triangularized  again  when  the  size  of  time  step  is 
changed. 

Although  this  procedure  is  tractable,  the  calculation  of  the  force  history 
vector  in  its  defined  form  is  computationally  unpractical  since  it  requires  the 
storage  of  all  past  incremental  displacements  and  the  formation  of  new  viscoelastic 
stiffness  matrices  at  each  time  step.  To  circumvent  this  difficulty,  a technique 
introduced  by  Taylor  (33)  can  be  exploited.  Taylor's  technique  requires  that  the 
relaxation  functions  in  Equation  B-3  be  in  the  form  of  a prony  series,  such  as 
Equation  B-2.  When  the  relaxation  functions  are  of  the  form  given  in  Equation  B-2, 
it  is  possible  to  restructure  the  force  history  vector  as  a summation  of  only  m^  + 
mQ  + 1 vectors  rather  than  a summation  of  all  history  vectors.  Moreover,  the 
exponential  form  of  the  relaxation  functions  permits  a simple  recursion  relationship 
for  updating  the  force  history  vector  at  each  time  step. 

Fluid  Foundation 

In  the  previous  section,  the  viscoelastic  boundary  value  problem  was  generated 
to  include  any  arbitrary,  axisymmetric  or  plane-strain,  geometrical  application. 
Now  the  geometry  is  specialized  to  plate-type  structures,  supported  by  a Winkler- 
type  fluid  foundation.  By  definition,  a Winkler  foundation  produces  a vertical 
resistance  force  proportional  to  the  vertical  displacement  of  the  nodes  at  the 
plate-fluid  interface.  Consequently,  these  additional  vertical  nodal  forces  at  the 
interface  can  be  expressed  as  [K^.]  {up, },  where  [K^]  is  a sparse  global  Winkler 
stiffness  matrix  with  non-zero  entries  associated  only  with  vertical  displacements 
at  the  plate-fluid  interface. 
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Considering  Equation  B-5  as  a force  equilibrium  equation,  the  influence  of  the 
fluid  foundation  can  be  added  directly  to  give: 

[K(Stn)]  (B-7) 

To  complete  the  inclusion  of  the  fluid  foundation,  the  components  of  the  Winkler 
stiffness  matrix,  [ K^l  , need  to  be  derived,  consistent  with  the  concept  of  virtual 
work. 

Superposition  of  Loads 

For  this  concept  the  geometry  is  assumed  to  be  an  axisymmetric,  half-space 
system  such  that  material  properties  vary  only  in  the  vertical  direction,  and  the 
radial  boundary  extends  to  infinity  (5).  Loading  is  represented  by  uniform  circular 
surface  pressure  discs,  and  their  location  is  completely  arbitrary.  Specified  load 
histories  are  assumed  proportional  to  one  another;  in  other  words,  the  pressure  of 
each  disc,  Pjg(t),  contains  a pressure  magnitude,  pj  , and  a loading  history,  g(t), 
which  is  common  to  all  discs.  With  the  above  assumptions,  the  state  of  stress  under 
any  pressure  disc  can  be  determined  by  solving  for  the  stresses  produced  by  that 
single  disc  and  using  stress  transformations  and  superposition  to  calculate  the 
stress  contributions  from  all  other  discs. 


USER  INSTRUCTIOM  FOR 
COMPUTER  CODE 

The  following  input  instructions  require  some  discussion  with  regard  to 
attempting  half-space  boundary  value  problems  and  using  the  mesh  generator.  By 
definition,  a half-space  boundary  value  problem  implies  that  the  boundaries  extend 
to  infinity  in  one  or  two  coordinate  directions.  However,  the  finite-element 
formulation  used  in  this  study  requires  selecting  boundaries  of  some  finite  extent. 
In  order  to  justify  a finite  boundary  representing  an  infinite  boundary,  it  is 
necessary  to  employ  St.  Venant's  principle  which  contends  that  if  the  finite 
boundaries  are  chosen  far  enough  away  from  the  area  of  interest,  then  the 
boundaries  and  boundary  conditions  will  have  a negligible  effect.  For  infinite  plates 
on  a fluid  foundation,  it  has  been  found  that  selecting  the  outer  radius  of  the  plate 
at  a distance  > 10  times  the  radius-of-relative-stiffness  is  sufficient.  For  a general 
half-space  problem,  several  attempts  of  extending  the  boundaries  may  be  required 
to  demonstrate  that  the  boundaries  are  not  influencing  the  solution.  The  general 
mesh  configuration  of  a plate  is  shown  in  Figure  B-1. 
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USER  INSTRUCTIONS  FOR  INPUT  AND  FORMATS  FOR  PROGRAM  VIS ICE 
BY  M.  G.  KATONA  (5) 


VIS ICE  IS  A FINITE  ELEMENT  VISCOELASTIC  PROGRAM  DESIGNED  TO  SOLVE 
FOR  STRESSES  IN  AN  AXISYMMETRIC,  VISCOELASTIC,  LAYERED  PLATE 
RESTING  ON  A WINKLER  FOUNDATION  SUBJECTED  TO  AN  ARBITRARY 
CONFIGURATION  OF  CIRCULAR  SURFACE  LOADS. 

A.  PROBLEM  INITIATION  AND  MASTER  CONTROL  CARDS 

CARD  Al.  {A6,11A6) 

COL  1 TO  6 MESHGN,  HOLLERITH  WORD  TO  CALL  MESH  GENERATOR. 

COL  7 TO  72  DESCRIPTIVE  TITLE  OF  PROGRAM. 


CARD  A2.  {2I5,3F10.0) 


COL  1 TO  5 
COL  6 TO  10 
COL  11  TO  20 
COL  21  TO  30 
COL  31  TO  40 


NUMAT,  NUMBER  OF  MATERIAL  LAYERS  IN  PLATE. 
NSTEPS,  TOTAL  NUMBER  OF  TIME-STEPS. 

DELT,  TIME-STEP  INCREMENT. 

ACELZ,  ACCELERATION  IN  Z-DIRECTION,  (GRAVITY). 
ANGFQ,  ANGULAR  VELOCITY  (ROTATING  PLATE) 


CARD  A3.  {4F10. 0,415) 


COL  1 TO  10 
COL  11  TO  20 
COL  21  TO  30 
COL  31  TO  40 
COL  41  TO  45 

COL  46  TO  50 

COL  51  TO  55 

COL  56  TO  60 


RLOAD,  LOADING  RADIUS  OF  WHEEL  (TIRE  PRINT) 

RMAX,  RADIUS  OF  PLATE 

PSI,  PRESSURE  OF  TIRE  PRINT 

GAMMA,  UNIT  WEIGHT  OF  SUPPORTING  FLUID. 

NELOAD,  NUMBER  OF  RADIAL  ELEMENTS  IN  CONTACT 
WITH  ONE  TIRE. 

NUMELR,  TOTAL  NUMBER  OF  RADIAL  SURFACE  ELEMENTS 
IN  THE  PLATE. 

NUMSPG,  TOTAL  NUMBER  OF  RADIAL  BOTTOM  ELEMENTS 
IN  CONTACT  WITH  FLUID.  (USUALLY,  NUMSPG-NUMELR) 
KBOUND,  OUTER  PERIPHERY  BOUNDARY  CODE... 


* KBOUND-0,  OUTER  EDGE  IS  FREE 

* KBOUND-1,  OUTER  EDGE  IS  HORIZONTALLY  FIXED. 

* KBOUND-2,  OUTER  EDGE  IS  VERTICALLY  FIXED. 

* KBOUND«3,  OUTER  EDGE  IS  COMPLETELY  FIXED. 

* 

* NOTE,  TO  REPRESENT  PLATES  OF  INFINITE  EXTENT 

* IT  IS  RECOMMENDED  TO  CHOOSE  KBOUND-3,  AND 

* RMAX-10*RADIUS  OF  RELATIVE  STIFFNESS 
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B.  MATERIAL-LAYER  CHARACTERIZATION  (ONE  GROUP  B FOR  EACH  LAYER, M) 


CARD  Bl.  (I5,2F10.0,2I5,2F10.0) 


COL  1 TO  5 
COL  6 TO  15 
COL  16  TO  25 
COL  26  TO  30 
COL  31  TO  35 
COL  36  TO  45 
COL  46  TO  55 


NELAY (M) , NUMBER  OF  ELEMENT  ROWS  IN  THIS  LAYER. 
HEIGHT (M),  THICKNESS  OF  THIS  LAYER. 

PO{M),  DENSITY  OF  THIS  LAYER.  (BODY  FORCES) 
IB(M),  NO.  OF  BULK  EXPONENTIAL  TERMS.  (4  MAX) 
IG(M),  NO.  OF  SHEAR  EXPONENTIAL  TERMS.  (4  MAX) 
BKO(M),  INFINITE-TIME  BULK  RELAXATION  MOD. 
GKO(M),  INFINITE-TIME  SHEAR  RELAXATION  MOD. 


CARD  B2.  (8F10.0) 


COL 

1 

TO 

10 

BK(MTYPE,1) , 

BULK 

COL 

11 

TO 

20 

BK(MTYPE,2) , 

BULK 

COL 

21 

TO 

30 

BK(MTYPE,3) , 

BULK 

COL 

31 

TO 

40 

BK(MTYPE,4) , 

BULK 

MOD.  OF  FIRST  E.XPONENTIAL . 
RELAX.  OF  FIRST  EXPONENTIAL. 
MOD.  OF  SECOND  EXPONENTIAL. 
RELAX.  OF  SECOND  EXPONENTIAL. 


CONTINUE  AS  REQUIRED  FOR  BULK  EXPONENTIAL  TERMS 


CARD 

B3. 

, (8F10. 

0) 

COL 

1 

TO 

10 

GK(MTYPE,1) , 

SHEAR 

MOD. 

COL 

11 

TO 

20 

GK(MTYPE,2) , 

SHEAR 

RELAX 

COL 

21 

TO 

30 

GK(MTYPE,3) , 

SHEAR 

MOD. 

COL 

31 

TO 

40 

GK(MTYPE,4) , 

SHEAR 

RELAX 

OF  FIRST  EXPONENTIAL. 

. OF  FIRST  EXPONENTIAL. 
OF  SECOND  EXPONENTIAL. 

. OF  SECOND  EXPONENTIAL. 


CONTINUE  AS  REQUIRED  FOR  SHEAR  EXPONENTIAL  TERMS 


C.  LOAD  FACTOR  CONTROL  CARDS  FOR  TIME  DEPENDENT  LOADING. 


CARD  Cl.  (I5,2F10.0) 


COL  1 TO  5 TLOAD,  LOAD-CASE  NO.  FOR  PROPORTIONAL  LOADING. 

COL  6 TO  15  Tl,  LOAD  DURATION  PARAMETER 

COL  16  TO  25  T2,  LOAD  DURATION  PARAMETER 


* TLOAD-1,  CONSTANT  LOAD  OF  DURATION  Tl. 

* TLOAD=2,  TRIANGULAR  LOAD,  PEAK  AT  Tl,STOP  AT  T2 

* TLOAD-3,  SINE-SQUARED  LOAD  PERIOD  OF  Tl. 

* 

* THE  ABOVE  LOAD-CASES  ARE  NORMALIZED  LOAD- 

* HISTORY  FUNCTIONS  WITH  A MAXIMUM  VALUE  OF  1 

* WHICH  MULTIPLY  ALL  EXTERNAL  LOADS  (PRESSURE 

* AND  POINT  FORCES).  FOR  EXAMPLE,  A CONSTANT 


liiASkldiiaiMiiyilHiMU 
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* LOAD  FOR  ALL  TIME  WOULD  BE  DESIGNATED  BY, 

* TLOAD=l,  AND  T1=DELT*NSTEPS . 

* 

* ADDITIONAL  LOAD-CASES  MAY  BE  PROGRAMMED  BY  THE 

* USER  IN  SUBROUTINE  TLOAD. 


D.  OUTPUT  CONTROL  CARD 


CARD  Dl.  (715) 


COL 

1 

TO 

5 

COL 

6 

TO 

10 

COL 

11 

TO 

15 

COL 

16 

TO 

20 

COL 

21 

TO 

25 

COL 

26 

TO 

30 

COL 

31 

TO 

35 

IPINT,  PRINT  INTERVAL  IN  TIME-STEPS. 

NDISPl,  FIRST  NODE  IN  A SEQUENCE  FOR  DISP. 
NDISP2,  LAST  NODE  IN  THE  SEQUENCE  FOR  DISP. 
NSTNl,  FIRST  ELEMENT  IN  A SEQUENCE  FOR  STRAIN. 
NSTN2,  LAST  ELEMENT  IN  THE  SEQUENCE  FOR  STRAIN. 
NSTRl,  FIRST  ELEMENT  IN  A SEQUENCE  FOR  STRESS. 
NSTR2,  LAST  ELEMENT  IN  THE  SEQUENCE  FOR  STRESS. 


* FOR  EXAMPLE,  IF  PRINTED  OUTPUT  FOR  EVERY  THIRD 

* TIME-STEP  IS  DESIRED,  SUCH  THAT,  DISPLACEMENTS 

* OF  NODES  10  THROUGH  20  ARE  PRINTED,  NO  STRAINS 

* ARE  PRINTED,  AND  STRESSES  OF  FIRST  5 ELEMENTS 

* ARE  PRINTED,  THEN  THE  ABOVE  PARAMETERS  WOULD 

* BE  INPUT  AS...  IPINT-3,  NDISPl-10,  NDISP2-20, 

* NSTNl-0,  NSTN2»0,  NSTRl»l,  NSTR2-5. 


E.  SUPERPOSITION  OF  STRESSES  CONTROL  CARDS 
CARD  El.  (15) 

COL  1 TO  5 INTSUP,  SUPERPOSITION  INTERVAL  IN  TIME-STEPS. 

FOR  EXA.MPLE,  IF  THE  SUPERIMPOSED  STRESSES  ARE 
DESIRED  OF  EVERY  FIFTH  TIME-STEP,  THEN  INTSUP-5 
INTSUP-0  RETURN  TO  MAIN  PROGRAM. 


CARD  E2.  (15,2F10.0) 

COL  1 TO  5 NLOADS,  NUMBER  OF  WHEELS 

COL  6 TO  15  PTOTAL,  TOTAL  WEIGHT  OF  AIRPLANE.  (FO.^CE) 

COL  16  TO  25  ZPLANE,  DEPTH  BELOW  THE  PLATE  SURFACE  WHERE 

SUPERPOSITIW  OF  STRESSES  IS  DESIRED  (NEGATIVE) . 
THE  ACTUAL  POINT,  (X,Y,Z)  IN  THE  PLATE  WHERE 
SUPERIMPOSED  STRESSES  ARE  CALCULATED  IS  THE 
POINT,  (0, 0, ZPLANE) . THE  VARIABLES  X AND  Y 
DESCRIBE  T-HE  LOCATION  OF  THE  WHEELS  ON  THE 
PLATE  SURFACE  WITH  RESPECT  TO  THIS  POINT  AND 
ARE  INPUT  IN  THE  FOLLOWING  CARDS. 

CARD  E3.  (3F10.0)  ONE  CARD  FOR  EACH  WHEEL  TO  BE  SUPERIMPOSED. 
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COL 

1 

TO 

10 

COL 

11 

TO 

20 

COL 

21 

TO 

30 

CARD 

E4. 

. (15) 

COL 

1 

TO 

5 

XLOC(N),  X-COORD.  OF  WHEEL  LOAD  CENTER 
YLOC(N),  Y-COORD.  OF  WHEEL  LOAD  CENTER 
PLD(N) , FRACTION  OF  TOTAL  LOAD  TAKEN  BY  WHEEL. 

EXIT  SUPERPOSITION  AND  RETURN  TO  MAIN  PROGRAM 

BLANK  CARD  SIGNALS  EXIT. 


* PRIOR  TO  CARD  E4.,  NEW  SUPERPOSITION  PROBLEMS 

* MAY  BE  DEFINED  AND  EXECUTED  BY  REPEATING  CARDS 

* El.,  E2.,  AND  E3.  THEREBY,  THE  USER  MAY 

* DETERMINE  THE  SUPERIMPOSED  STATE  OF  STRESS  AT 

* NUMEROUS  POINTS  AND/OR  FOR  DIFFERENT  WHEEL 

* CONFIGURATIONS  AS  LONG  AS  THE  BASIC  TIRE-PRINT 

* RADIUS  IS  THE  SAME. 


F.  TERMINATION  CARD  (THIS  CARD  FOLLOWS  LAST  DATA  DECK) 


CARD  FI.  (A6) 

COL  1 TO  6 STOP  HOLLERITH  WORD  TO  STOP  EXECUTION. 
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RMAX 


Figure  8-1.  Mesh  configuration  of  a piale  structure. 
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